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Recent Results in Signal Processing and Communications 


Foreword 

Phis special issue consists of ten papers which have been selected from among thosi 
jresented at the Conference on Signal Processing, Communications and Networking, hel( 
It the Indian Institute of Science, Bangalore during July 16-19,1997. These papers repor 
:ome new results in the fields of signal processing and communication. They cover a widi 
ipectmm of topics of current interest, namely, parameter estimation in the presence o 
R.icean fading, transformations for mapping circular array manifold to that of unifom 
inear array and the associated errors, use of fractal models and multiresolution frameworl 
n image reconstruction and interpretation, group codes with rotational invariance property 
ise of group delay function in power spectrum estimation of complex signals, analysis o 
stochastic gradient lattice algorithm, adaptive notch filters, schemes for I-Q imbalanci 
jorrection, application of cosine-modulated filter banks in reconstruction of periodic erro 
jursts in oversampled data and performance analysis of two versions of TCP over mobiL 
•adio links. 

Use of sensor arrays in the field of mobile communication to improve the performanci 
rf a cellular system has been an active area of research. Recently, maximum likelihoo( 
md least squares based techniques have been proposed for estimating direction-of-arriva 
POA) of the received signal from a mobile unit and the associated angular spread due ti 
icattering in the presence of Rayleigh fading using a linear array. The first paper by Hai 
md Otterson addresses a similar problem for a Ricean fading channel. 

Recently, several high resolution DOA estimation methods and the effect of certaii 
jre-processing schemes, known as smoothing techniques, on the performance of thes- 
nethods have been studied by many authors. Most of these studies have been restricte( 
;o uniform linear arrays because several of the high resolution methods and smoothini 
schemes are applicable only to such arrays. However, in practical applications where 360 
coverage is required, we use uniform circular arrays. One can apply a transformatior 
Dased on phase mode excitation concept, to the data from the circular anay and use th 
techniques applicable to uniform linear arrays. The second paper by Maheswara Redd 
md Umapathi Reddy discusses the errors arising out of this transformation and their effec 
an the performance of the DOA estimation. 

In computer-aided tomography, the objective is to reconstruct a cross-section of an objec 
From measurements that are strip integrals of some property of the object. A popula 
technique for image reconstruction is the filtered back projection. The third paper b 
Chowdhury, Barman and Ramakrishnan addresses the problem of reconstruction fror 
floisy projections, restricting their discussion to parallel beam tomography. One of th 
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advantages of their approach is the reduction in the time of exposure to the p 
radiation. 

Digital communication over additive white Gaussian noise channel can be mo< 
transmission of a point from a finite set of points of a finite-dimensional vectc 
The problem of signal set design then reduces to choosing a specified number c 
in a space of specified dimension in such a way that the minimum distance bet\ 
points is maximum. One can obtain good signal sets in large dimension starting 
signal set of small dimension using group codes. The fourth paper by Jyoti Bali an( 
Rajan reports rotational invariance properties of coded signal sets obtained by t 
constraction for a class of two- and four-dimensional signal sets. 

Spectrum estimation is a highly researched field. The basic periodogram spec 
mate (PSE), which has good resolution, low bias and gives good signal detectabilib 
from large variance. A group delay approach, with modification for the group de 
serves the advantages of PSE and reduces the variance of the estimate. The fif 
by Narasimhan, Plotkin and Swamy applies the group delay function with mod 
to complex signals to estimate the power spectrum, and discusses its applicatioi 
Wigner-Ville distribution. 

Adaptive HR notch filters provide superior performance at lower computatic 
relative to their FIR counterparts in suppressing narrowband interference and p 
line enhancement. A diverse choice of biquad designs are proposed for perforr 
above functions by several authors. The sixth paper by Krishna and Hiremath pr 
lucid presentation on how all these designs are equivalent in their asymptotic perfc 

The adaptive lattice filter has several desirable characteristics such as orthogon; 
of the input, faster convergence and simple stability checking criteria. The st 
gradient lattice algorithm is commonly used for adapting the FIR lattice filter coe 
because of its low computational complexity. The seventh paper by Negi and P 
presents some theoretical results on bias in reflection coefficient and convergem 
algorithm. 

In a coherent receiver. In-phase and Quadrature phase signals are derived by dei 
ing the IF signal. Practical demodulators give rise to some amplitude and phase im' 
in I-Q signals. The eighth paper by Sivannarayana and Veerabhadra Rao suggests 
in time and frequency domains for I-Q imbalance correction. 

The restoration of erasure bursts in oversampled multimedia data has assumed ii 
significance in the context of cell loss in ATM networks. Many heuristic method: 
linear interpolation, repetition, muting, and systematic approaches have been co 
in recent years. The ninth paper by Jayasimha and Hiremath presents a low-coi 
approach, based on cosine modulated filter banks, for reconstructing periodic 
bursts. 


Foreword 
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versions of TCP (the original, and the Tahoe version) over a Rayleigh fading radio link; the 
analysis reveals the effect of vehicle speed, and the average SNR required for achieving 
adequate throughputs. 

We would like to express our thanks to Prof. N. Viswanadham, Editor, Sadhana, foi 
inviting us to bring out this special issue, and also special thanks to the editorial staff o1 
the journal for their effort and cooperation in bringing out this issue. 


February 1998 


V UMAPATHI REDD\ 
ANURAG KUMAF 
Guest Editor; 



)ddhana, Vol. 23, Part 1, February 1998, pp. 5-15. © Printed in India. 


Parameter estimation using a sensor array in a Ricean 
fading channel 

K V S HARI' and BJORN OTTERSTEN ^ 

^Department of Electrical Communication Engineering, Indian Institute oi 
Science, Bangalore 560 012, India 

^Department of Signals, Sensors and Systems, Royal Institute of Technology 

S-100 44 Stockholm, Sweden 

e-mail: hari@ece.iisc.emet.in; otterste@s3.kth.se 

Abstract. The estimation of the Direction-Of-Arrival (DOA) and the variance 
of the angular spread, using an array of sensors in the case of a Ricean channei 
is considered, using the Maximum-Likelihood, Least-Squares and Weightec 
Least Squares criteria. The Cramer-Rao bound is also obtained for the problen 
of interest. Simplification of the cost functions to reduce the dimension of the 
problem has been carried out and the performance of the methods has beer 
studied based on numerical experiments. 

Keywords. Antenna arrays; fading channels; spatially distributed sources. 


1. Introduction 

Fhe use of sensor arrays in the field of mobile communication to improve the perfor 
mance of a cellular system has been an active area of research recently. An importan 
issue is to efficiently use the available channel bandwidth to provide services to as manj 
users as possible. Techniques using an array of sensors have been proposed which esti 
[uate the Direction-Of-Arrival (DOA) of the received signal from a mobile unit and the 
associated angular spread due to scattering, to form multiple beams on the same chan 
ael and increase user capacity (Yeh & Reudink 1982; Anderson et al 1991; Balaban & 
Salz 1992; Ohgane et al 1993; Zetterberg Sc Ottersten 1995). Estimation of DOA anc 
angular spread of scattered field using a Uniform Linear Array (ULA) has been carriec 
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criteria. A comparison of the performance of these methods is carried 
simulations. 


2. Data model 

In a mobile communication scenario, the narrowband signal from 
received at the base station using a sensor-array, is assumed to be 
large number of signals with different strengths and arriving from an| 
to the direction of the source. This model has been verified by experii 
the scattering effects of the channel (Adachi et al 1986). In this report 
with an inter-sensor distance, dx (in wavelengths) is considered. The: 
to a unit-amplitude narrowband signal impinging from a direction 9 
broadside of the array), is known as the array response vector for 
denoted by a(0) and its I:th element is defined as' 

[a(0)]A: = expljlrrdxik - 1) sin(6i)]. 

Assume that a single source is transmitting a narrowband signal s 
9 which is being received by the ULA due to scatterers in the vicinil 
noiseless output of the fcth sensor as a function of time t, can be writ 

ykit) = 

where gn(t), otnit) are the amplitude and phase factors due to the ni 
the deviation with respect to 9 due to the nth scatterer. Stacking the o 
into a vector, y(0, the snapshot of the array output can be written as 

y(r)=h(r)r(0, 

n 

where h(f) denotes the channel response vector at time instant t. 

2.1 Fading channels 

A channel is said to be & fading channel if the amplitudes and phasi 
scatterers and the directions of the scatterers are random and vary with 
Braun & Dersch 1991). Further, if these are independent and identica 
the number of scatterers is very large, the Central-Limit Theorem cai 
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ignal amplitudes, phases and the directions needs to be done. Trump & Ottersten (199€ 
lave assumed that 

E[gn{t)e^‘'''^'^gm{s)e~j°‘"'^^h = Q, n^m, t f^s, 

= 1, Otherwise, 

£[e^'“”W]=0. 


iased on these assumptions, Cou[h(t), h(T)] = Ch8(t — r). Many distributions c 
have been proposed by Anderson et al (1991), Parsons & Turknani (1991), Proaki 
1991) and Trump & Ottersten (1996) and in this paper it is assumed that is small an 
)„ ^ Af(mg, ag). 

1.1a Rayleigh fading channel: In an urban environment, usually, there is no direc 
lath between the mobile unit and the base station sensor array and it is assumed th£ 
n/i(r) = £'[h(r)] = 0 (Zetterberg & Ottersten 1995; Trump & Ottersten 1996) and sue 
i channel is known as a Rayleigh fading channel (Proakis 1989). Assuming, mg =0,C 
s given as (Trump & Ottersten 1996) 

Q = Rfl • B(6>, ag), 

vhere the kith element of B is given by 

[B]« = exp{-2[7r(^ - l)dxfag cos^(6i)}. 


ind 


R« = a(0)a^(0). 

\4odel Ml - Let s(t) be a random signal uncorrelated with h(0 and ^’(0 = a exp [7 (ps (t] 
vith a being a deterministic quantity and (psif) being distributed uniformly between 0 an 
ijr. Then m3;(r) = 0, Cy{t) = Ch\cc\^^ 

It can be easily shown using elementary probability theory, that y(r) ~ ^/'(O, Cy{t] 
br model Ml. 

Remark 1, It is to be noted that if s{t) is Gaussian, y{t) is no longer Gaussian. 

Model M2 - If s{t) is a deterministic signal, then my(t) — 0, Cy(t) = Ch |s(^)p an 
m - ^( 0 , Cy(t)), 


12 Ricean fading channel 
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Case 1: Let y be a deterministic but unknown quantity. 

Case 2: Let y ~ iV (yo, ay) independent of other parameters. Since i 
assumed, yo # 0. It is clear that the new channel response vector in be 
and has a non-zero mean of the form ya(0). This channel is said to 
channel because the envelope has a Rice distribution (Proakis 1989). 

Model M3 - Let s{t) be a random signal uncorrelated with h(t) anc 
Ml. Then m^Ct) = 0, Cy{t) = (Q -f mhmf) lap.As in the case o 
Gaussian. 

Model M4 - s{t) is a deterministic signal. Then m^, (t) = m;, 5(r), Cy 
y(r) ~ Afiaiyit), Cy(t)). 

2.3 Noisy data 

Consider the noisy array output as 
x(t) = y(t) + nit), 

where n(f) is the additive noise vector which is zero-mean Gat 
= 0, £[n(f)n'^(r)] = a„I and £:[n(f)n'^(T)] = 0. Then C 
for any of the above models. In this paper, the Ricean fading channel 
M3 and M4) with y being deterministic but unknown. The covaria 
matrices for x(t) belonging to these models can be obtained as 

= CTy (Js Ra -+■ CTj Ri) -j- cr,;l, 

Cx = tTfRi -H o>iI, 

where Oy = |y |^, Os = and mjc is the mean of x(0. 

3. Parameter estimation 

It is clear that there are two cases of interest: (i) Random source (i 
deterministic source (model M4). The estimation of parameters for 
channel assuming a random source (model Ml) was carried out (Trum] 
In this paper, the estimation of the parameters for models M3 and M4 
Likelihood and Weighted Least Squares criteria is presented. 

Problem statement. Given N snapshots of the array outputs, x(/i), x(t2 
the parameter vector n = \6, ao, cr?, cr„, Ov]^ for model M3 and • 
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( 


vhere is the Fisher Information Matrix. For model M3, the kith, element of is (Trumj 

k Ottersten 1996) 


vhere d{.)/drjj^ denotes differentiation with respect to the kth parameter of rj. For mode 
AA, the kith element of Jjj can be easily shown to be 


[J^]w = ^Tr 


^- 1^-1 2C~^ 


vhere C;c and R;c denote the mean, covariance matrix and the correlation matrix o 
tit). 


>. Maximum-likelihood estimation 

jiven the Gaussian nature of x(ri), xfti), - • •, x(tAr), the negative log-likelihood functioi 
or model M3 is given as 

iMliO, OQ, Os, On, Oy) = log(det(R;c)) + Tr[R“^R] 
ind the conditional negative log-likelihood function for model M4 is given as 

iMii^, oe,o„,Oy,sit), t = l,...,N) 

= log(det(C;c))4-Tr[C;’ M], 

vhere M is defined as 

M = R -1- mm^ — 2mm^ 

vith 

1 ^ 

R = -/£xit)x^it), 

^ tl 

1 ^ 

m = — Yxit). 

R is the data correlation matrix^ and ih is the sample mean of x(f). The maximun 
likelihood (ML) estimate of rj is obtained by minimizing 1ml in the parameter space of i 
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6. Eeast squares estimation 

As the solution to the ML problem is computationally expensive, a look at ( 
criteria like the least squares criterion is worthwhile. 

6.1 Weighted least squares 

The general form of least squares is the weighted least squares (WLS) cost fu 
can be expressed as 

I = Imjsiv) = I|WH/2(r^ _ R)w'/2||2 

= Tr[(R;c-R)W(R;,-R)W], 

where W is a positive definite weighting matrix. The choice of W is usuall 
that the error-covariance of the parameter vector, 77, is minimized. Denoting 
of 770 77, the error in the parameter vector is given by 

-1 

fi = (vo-‘n) ^ -H '(770) 

where H is the Hessian of the cost function given by 

dH 

[H],7 = 

drtidrjj 

The fth element of the gradient vector /, of the cost function can be obtainec 

r = ^ = 2Tr[(R;c-R)D,], 

drji 

D,=W^W. 

dm 

Following the development by Trump & Ottersten (1996), 

As Af(0, Q), 

where 

Q= lim NEll'il'fl 

N-^oo 

Hence, the asymptotic distribution of the estimation error is given by 
ViVi) ~ As ATfO, C), 
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rhe yth element of Q is given as 

= 4£[Tr{(R;,-R)D/} 

xTr{(R;,-R)D,}] 

= 4Tr{R;cD,}Tr{R;,D,} 

-4Tr{R;cD/}£[Tr{R^D,}] 

-4£[Tr{RDi}]Tr{R;,D,} 

+ 4£:[Tr{RD,}Tr{RDy}] 

= 4£[Tr{RDi}Tr{RDy}] 

-4Tr{R;cD,}Tr{R;,D^} 

as Em = R. The second term in the above equation can be written as 
= 4£[Tr{RD/}Tr{toj}] 

= E E[xf(t)Xrn{t)xlit)Xp{t)mi]op[Dj][„. (2 

txlmop 


dl dl 

drji drij 


6.1a Random source (model M3): For the random source case, since nijcCO = 0, th 
above product of four Gaussian random variables with zero mean can be expressed as 

= E {E[xUt)^mit)]E[xl{x)Xp{x)] 

txlmop 

+ £:[X*(r)Xp(T)] E[xl(x)Xm(t)])[Di]op[Dj]lm 
= 4Tr{R;cD,}Tr{R;cDy} + ^Tr{R^DmBj}. 

Thus 


dl dl 


drii dr\ 


'JJ 




3R;c, 


= -Tr R^W^WR;cW^W 


N 


9?7i 


dr)j 


It was shown (Goransson 1995) that W = R^ * would yield 


[C-% = Tr 


9 Ra: p_i 9 R;c p_i 

-TIT 


drjj 


dr)i 
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Differentiating w.r.t. cry, setting to zero, snd simplifying, tli 

hold, 

an<Ty(Ts + r3tl2ds + = b\, 

a2i(yyO's +a22<^s + aiicr„ =^2, 
a220'y CTy + a320's + «12Crn = 

where an = Tr[R~*RaR’'M^, ai2 = Tr[R 'R^jR 
Tr[R“^RflR'“'Ra],a22 = Tr[R"’R;,R“‘Rf, 1,^32 = Tr[R"'R/,R 
A2 = Tr[R-^Ra], h = Tr[R-’R6], Using the above equations, 0 
sions for the three parameters as 

- ^22^1 ~ '^11^3 ~ g>[fll2r^22 - ^32<3n ] 

" [< 313<222 — < 2 l 2 < 2 ll] 

„ _ an^i — Qi3^2 ~ (tslanan - ^ 13 ^ 22 ! 

^ crjaf]-ai3a2i] 

^ Cl + C2 + C3 

0's — ? 

C4 


where 


ci=hi 1- 




[( 2 ^, - ai 3 ( 32 l] [' 5 > 13 « 22 -<;il 2 «ill/ ' 


C2 = b2 


flisaii 


[afj -ai3<32l] 


. C3 = h 


«I 3«11 


[< 213^^22 ~ <^l\ 2 <^\\ J 


( n'il[aiiai2 — ai3<222] ^I3[ai2<222 “ ^3'>n 

an - -2 -^- 7 -^ 

Kl - « 13 « 2 l] f<^13«22 - «I2"I I 

The weighted least-squares cost function can now be recast as 


^WLs(^, <ye) = Tr(((Ty(7yRa + (7,,R/, -I- a„I)R ' - I)“], 


and the search is now over a two-dimensional space of [0, ao] whic 
less expensive than before. 


6. lb Deterministic source (Model M4): For the model M4, the ra 
have non-zero means and thus 


Parameter estimation using a sensor array in a Ricean fading channel 


1 


\fter some tedious calculations, the ijih element of Q for this case is given as 

[Q],7 = ^ Tr{R;, D,} Tr{R;, + V D/} Tr{(R;)T D,} 

-8 Tr[m;, m« D,] Tr[m;c mH D;], 

where 

r; = £[x(0xT( 0.] = m;,(0mj(0. 

Remark 2. 

9 It is difficult to obtain W from the above equation for Q, which minimizes C. 

9 The same result holds good for model M2 also, with the appropriate correlation matri> 


5.1c Least squares: A more popular criterion, which is simpler, is the least square 
criterion, defined for W = I as 

/ls = Tr[(R - R^)(R - 


Random source (Model M3 J - As in the case of the WLS criterion, since the cost function i 
quadratic in as,ay,(7n, these parameters can be separated. Using some simple identities, 
one can obtain 


CTn 


;6Tr[R] - LTrlRRb] - asL(p- Tr[RfcRH]) 
L{p - L) 


CTy 


Tr[RRJ-Tr[R]^a,(^~L) 

crsL(L - 1 ) 

numerator 




denominator ’ 


where 


numerator = Tr[R]((L2 - 2L + - L){L - 1))) 

-Tr[RRJ(l/(L - 1)) +Tr[RRi](L/(y3 - L)), 
denominator = L - (L(/8 - Tr[RfeRt]))/(;6 - L) - - L)/{L - 1). 


The least-squares cost function can now be recast as 
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Table 1. Mean squared error (deg^) in the DOA (=10 deg) vs L for erg — 1 


L 

ML 

LS 

4 

1.2560e-02 

1.2512^-02 

6 

4.4170e-03 

4.3335e-03 

8 

3.7610i?-03 

3.7726e-03 

10 

2.3106^-03 

2.2610e-03 


7. Numerical study 

An experiment to study the performance of the algorithms based oi 
criteria, is presented next. 

Experiment. A scenario with 9 = 10'’, Os = 10, cr„ = l,cry 
ULA is considered. Various values of ae,L are considered as i 
L = 4,6, 8, 10,12. 

The estimates of ^ and ae are obtained for each of the above comb 
the ML, LS and WLS criteria using the same data vectors. In this r 
search method is used to obtain the parameters. The updated vec 
given by 

fj(k 4-1) = ^{k) - ix(k) g, 

where k denotes the iteration, H is the Hessian and g is the gradie 
considered and p(k) is the step-size at the feth iteration. The initial e 
by using the ESPRIT algorithm and the initial estimate of ag = 0 
sample statistics of the estimates are obtained from 200 independe 
Effect of number of sensors: Table 1 presents the MSE in the DO.;* 
a particular value of the angular spread, ag. 

• The MSE decreases as L increases for all methods which agrees 

• For any value of L, the performance of the LS method is very ( 
method while the WLS method performs poorly in comparison v 
This could be due to use of the estimate of the correlation ma 
weighting matrix. 

Table 2 presents the MSE in the DOA as a function of the angular 
for a particular value of L. 


Table 2. Mean squared error(dee^) in the DOA (=10 deg) vs cta for L = 10 
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i It is clear that the performance deteriorates as the angular spread increases, for al 
methods. 

» As observed before, WLS performs poorer than LS and ML while the performance o: 
ML is the best among the methods. 

S. Conclusions 

Estimation of parameters, the DOA and the variance of the angular spread, using an arra] 
3f sensors in the case of a Ricean channel is considered, using the maximum-likelihood 
[east-squares and weighted least squares criteria. The Cramer-Rao bound is also obtainec 
for the problem of interest. Due to the quadratic nature of the least-squares criteria, sim 
plification of the cost functions to reduce the dimension of the problem has been carriec 
Dut. The performance of the methods (in terms of the mean-squared error in the estimate: 
of the parameters) has been studied based on numerical experiments which show that th( 
maximum-likelihood and least-squares methods perform comparably while the weighted 
[east squares method is slightly poorer than the other methods. This could be due to th( 
use of an estimated correlation matrix as the weighting matrix instead of the true one. 
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Spatial smoothing with uniform circular arrays 


K MAHESWARA REDDY * and V U REDDY 

Electrical Communication Engineering Department, Indian Institute of 
Science, Bangalore, India 

* Present address: Centre for Aeronautical Systems Studies and Analysis, De¬ 
fence Research and Development Organisation, New Tippasandra, Bangalore 
560075, India 

e-mail: mahesh@cassa.emet.in; vur@ece.iisc.emet.in 

Abstract. In this paper, we extend and analyse spatial smoothing with uni¬ 
form circular arrays (UCA’s). In particular, we study the performance of the 
Root-MUSIC with smoothing in the presence of correlated sources, finite data 
perturbations and errors in transformed steering vector that arise due to some 
approximations made to enable the extension of the Root-MUSIC and smooth¬ 
ing to UCA. Expressions are derived for the asymptotic performance of the 
Root-MUSIC with smoothing applied to the transformed UCA data. An attempt 
has been made to bring out the impact of both forward and forward-backward 
smoothing. Computer simulations are provided to demonstrate the usefulness 
of the analysis. 

Keywords. Uniform circular arrays; phase modes; direction of arrival esti¬ 
mation; spatial smoothing. 


1. Introduction 

Uniform circular arrays (UCA’s) are commonly employed when 360° coverage is required 
in the plane of the array. Circular arrays are non-uniform linear arrays, and hence, the 
rooting techniques and preprocessing schemes like spatial smoothing cannot be directly 
applied to these arrays. Tewfik & Hong (1992) showed that it is possible to extend the 
Root-MUSIC to UCA using the phase mode excitation concept. Mathews & Zoltowski 
(1994) proposed real beamspace MUSIC to UCA (UCA-RB-MUSIC) which yields re¬ 
duced computation and better resolution. They also studied the direction of arrival (DOA^ 
estimation performance of the UCA-RB-MUSIC. 

While extending the rooting techniques to UCA, all the authors assumed that some of the 
terms in the transformed steering vector of UCA, where the transformation is performec 
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DOA estimates obtained with the Root-MUSIC even when the number of snap 
to infinity, and we analyse the effect of these errors in this paper. We also exi 
smoothing to UCA and analyse the effect of smoothing in the presence of correla 
and finite data perturbations. We discuss the impact of both forward and forwarc 
smoothing. 

In § 2, we provide a brief background. We propose in § 3 forward spatial 
highlighting the assumptions made for extending the Root-MUSIC and spatial 
to UCA and the errors associated with these assumptions. In § 4, we analyze 
mance of the Root-MUSIC with forward and forward-backward smoothing ap 
transformed UCA data. We present the results of computer simulations in § 6 ar 
the paper in § 7. 

2. Background 

Consider a UCA with L identical and omni-directional sensors. Let r be the n 
array and d be the circumferential spacing between the elements. Let 9 denot 
(azimuth angle) measured in the plane containing the elements. We assume foi 
that the sources are in the same plane as the UCA. The steering vector of the 
the centre of the array can then be expressed as ’ 

2^^{9) — , ^i?cos(6l-2:nr(L-l)/i:)jT 

where § = 2nr/X, X is the wavelength and (.)^ represents the transpose of (.) 

Consider the phase mode excitation of the UCA. The weight vector that 
array with mth phase mode is given by (Davies 1983) w,^ = 

^j2:7rm(L-i)/Lj •pjjg ajj-^y pattern for the mth phase mode can be shown to I 
(Davies 1983; Mathews & Zoltowski 1994) 

/m(0) = w«a,(0) = JH(^)e>'' 

OO 

+ -V<m . 

q=\ 

where V is the maximum number of phase modes given by (Davies 1983) V : 
Jm{^) is the Bessel function of the first kind of order m, h = Lq + m, g = Lc 
represents the complex conjugate transpose of (.) and Ixj denotes the largest 
than or equal to x. The first term in (2), the principal term, becomes dominant 
than 0.5A. In our analysis, we consider d < 0.5X and assume the second term 
negligible. 

The normalised transformation matrix F to excite the array patterns corres 
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vhere 

= x/Idiag[7x.(f). J^m, 

a(e) = ^-^-(^-1)0,..., 1,..., ^jT>6f^ ^4, 

ind Aa(0) is the contribution due to the second term in (2). Note that the vector a(0) 
las a structure similar to that of the steering vector of a uniform linear array (ULA). This 
suggests that we can extend the spatial smoothing to UCA provided the term Aa(i9) is 
legligible. We treat Aa(0) as the error in the transformed steering vector, caused due tc 
ipproximation. 

J. Forward spatial smoothing with UCA 

A.ssume that M sources are impinging on the UCA and the DOA’s of these sources are 
9i, 02 , ■■■, Om- If we assume that the signal and noise are uncorrelated and noise is 
spatially white with variance then the covariance matrix at the output of UCA can bf 
jxpressed as 

Rc = AcSA” + ah, (5; 

where S is the signal covariance matrix, I is an identity matrix and Ac is the matrix o 
direction vectors of the UCA. From (5) and (3), the covariance matrix that we obtain afte: 
applying the transformation F can be shown to be 

R' = F”R,.F = Jf A S A”J| + ah + AR, (6 

where A = [a(0i),..., a(0A/)] and AR = AA S A^J^ + J|A S AA^ + AA S AA^ witl 
AA = [Aa(0i),..., Aa(0M)]. Note that the size of R' is (I'D +1) x (2P+1), and hence 
the number of sources (M) should be less than (2D + 1) for any MUSIC-type algorithn 
to be applied to R'. 

Spatial smoothing is a preprocessing scheme originally proposed for ULA to alleviati 
the ill effects of correlation. This scheme can be extended to UCA by applying the transfer 
mation to the covariance matrix R^ provided the term AR is negligible. If the source 
are in the same plane as the UCA or if all the sources are at the same but known elevatioi 
angle, then is a known matrix. 

Let K be the number of virtual subarrays (since the subarrays are not physically avail 
able). Then, the forward smoothed covariance matrix R^ is given by 
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wheretheprewhiteningmatrixRnu, — [^ Z^J| Z/] 

equal to the number of coherent signals present, then the forward sp 
up the rank of the smoothed signal covariance matrix and makes it r 

1985 ). 

Using the structure of a(0) as given in (4), it can be shown (Red 
effective correlation coefficient (p/) between the sources after fon 
case of UCA, is given by (assuming AR to be negligible) 


\Pf\ = 


s.in{K{ei - Bj)IT) 

^ K sin((0i - dj)l2) 


where d-, and 0j are the DOA’s of the i th and y th sources respectively, 
coefficient between these sources before smoothing. Note from (9) tl 
on the individual directions of the sources, but is dependent only on tl 
If this angular separation is 90°, then pf becomes zero for AT=4. O 
separation is 180°, then two subarrays (if=2) are enough to force 
from (9) that pj is independent of the spacing between the elemer 
less than 1/2 (making Aa(0), and hence, AR to be negligible), 
forward spatial smoothing is applied to ULA, the effective correlatic 
is dependent on the individual directions of the sources and also o 
the elements. 


4. Performance of the Root-MUSIC with smoothing 

In this section, we analyze the performance of the Root-MUSIC witl 
backward smoothing applied to the transformed UCA data. 

4.1 Forward spatial smoothing 

Consider the smoothed covariance matrix after prewhitening (sec 
with (6) and (7), we obtain 
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md 

1 ^ 

^ /=i 
1 ^ 

= + ASAA«Jj-']Z,rH^. (12 

5/ is the smoothed signal covariance matrix and A/ is the virtual subarray direction matrix 
[n writing the RHS of (12), we assumed A A to be small and neglected the terms containinj 
nore than one AA. Note that Rf is the smoothed covariance matrix that we would get i 
Aa(0) (see (3)) is zero. In practice, however, this term may be small but non-zero, thereb; 
resulting in errors in the DOA estimates when we apply the Root-MUSIC to (Ry)u;. Wi 
low analyse the effect of this term (i.e., ARp) and that due to finite data perturbations oi 
;he DOA estimates. 

Let (Rj-)w denote the estimated covariance matrix from finite number of snapshots. Thi 
:an be expressed as 

(R})u, = (R})u, + ARp = Rf + ARf + ARp, (13 

where ARp represents the perturbation due to finite data. Note that ARp is random whil 
ARf is deterministic. If we assume that the noise at the output of the sensors is comple: 
circularly Gaussian distributed, then the mean square error (MSB) in ith DOA estimate 
lue to both the finite data perturbations and the error due to approximation (i.e. due t 
ARf), can be shown to be (Rao & Hari 1990) 

ElAOh — + 2[Rg(a^ARF)3)]^ 

'■ 2[vH(0,)RH^P„R„^v/,(0,)]2 

1 ^ ^ 

E E 

p=lq=\ 

K K 

= (15 

p=\q=\ 

a = PnRuvvVyi (0;); f = (RF)^RnHiy/'(0i)! 

(Rf). = R„v.A/S/A}^R«^, (16 

where = E[yp(r)y^(r)], yp(r) is the output vector obtained from the pi 

virtual subarray after prewhitening, N is the number of snapshots, Fn is the projectio: 
matrix onto the noise subspace of Rf, v/, (9) is the derivative of Vf (0) w.r.t. 9 with v/(6 
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Reddy & Reddy (1996a) showed that the smoothing reduces the nois 
to finite data in addition to reducing the correlation among the imping 
As the number of snapshots tends to infinity, the MSB in the DOA 
because of the error due to the approximation (cf. (3)). This error, whi( 
asymptotic error in the DOA estimate, is deterministic and given by 

Note that this error increases as d tends to X/l since AA becomes lai 
values of d. Let us first assume that the sources are uncorrelated. Then 1 
can be shown to be (see Reddy & Reddy 1996b) 


1 [j?g(aHR„v,Aa/(g,))]^ 

^o[vH(0/)RH,P„R„wV/,(0,)P’ 


where 

Aa/(0,-) = i 

^ /=l 

which we define as the effective error along the direction of the i th sourc 
steering vector due to the approximation. Note from (19) that the asy 
ith DOA estimate is dependent only on the effective error vector Aa/(6 
vector can be shown to decrease with spatial smoothing (see Appendix ^ 
1996b). Thus, we can expect the smoothing to improve the asymptotic 
Root-MUSIC applied to the transformed UCA data. Expression (19), 1 
only when the sources are uncorrelated. For the case of correlated sourc 
(18) leads to lengthy expressions even for a two-source case, and henc 
here. 

4.2 Forward-backward spatial smoothing 

Forward-backward spatial smoothing (William et al 1988) can also b 
by applying the transformation ' to the covariance matrix R^ (given 
term AR is negligible. If K is the number of virtual subarrays, then the 
smoothed covariance matrix is given by 


shown (Williams et al 1988; Pillai & Kwon 1989) that the forward-backward smoothing 
;an handle up to [2(2D -I- 1)/3J coherent signals in contrast to the forward smoothing 
ivhich handles up to [(I'D -I- 1)/2J coherent signals only. 

Taking into account the structure of a(0) given in (4), it can be shown that the effective 
correlation coefficient (pfb) between the sources after forward-backward smoothing, ir 
Jhe case of UCA, is given by (assuming AR to be negligible) 


Ififbl = \P\ 


sin (/if {9i — 9j)/2) cos ijr 
K sin((0; - Oj)/!) 


( 22 ; 


where p — \p\e^'^. Observe from (22) and (9) that the effective correlation with forward- 
backward smoothing is same as that with forward smoothing only, when i/^ is zero, i.e., f. 
is real. When ijr is an odd multiple of nil, the effective correlation with forward-backwarc 
smoothing reduces to zero for all K. 

The mean square error (MSB) in the ith DOA estimate that we obtain by applying th( 
Root-MUSIC with forward-backward smoothing to the transformed UCA data, due to bott 
the finite data perturbations and the error due to approximation (cf. (3)), can be shown tc 
be (see Rao & Hari 1990) 

“ iV/i:22[vH(0,)R«.P«RwV/, 

K K 
p=lq=l 

+ 2[/?e(a«ARFBjS)]^]. (23 

0! = P^Rnwy/i (^i); P — (Rfb)^ Rnwy/'(^i)i (24 

('ll 'I n \ • Ml — I(ARF)*i 

(RFB)j = RnwA/lS/i!,A^K„^^, AKfb = -^, (2b 

where and are the noise projection matrix and the signal covariance matrix, respec 
tively, obtained with the forward-backward smoothing. The first two terms in (23) are du( 
to finite data perturbations and the third term is because of the error due to the approxima 
don. The performance with forward-backward smoothing is the same as that with forwar< 
smoothing when S is real. Hence, the expression (19) can be used for the asymptotic per 
formance of the Root-MUSIC with forward-backward smoothing also, when the source 
are uncorrelated. 

The analysis carried out, so far, assumes omni directional sensors. Spatial smoothing 
can also be extended to UCA with directional sensors (see Reddy & Reddy 1998) and thi 
analysis of Section 4 is applicable for UCA with directional elements. 


5. Numerical and simulation results 
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Figure 1. Norm of the effective error vector in the transformed steering vec 
a function of circumferential spacint (L = 30). 


vector was the sum of noise and signal vectors which were generated sep 
vector consisted of zero mean, unit variance, independent complex cii 
random variables. In the case of uncorrelated signals, the signal vector 
zero mean, independent complex circularly Gaussian random variables 
af, where a]' was chosen to give the desired signal powers. In the case 
nals, the second signal S 2 (t) was generated as S 2 it) = psi (t) + y/l — \p 
and s{t) are zero mean, independent complex circularly Gaussian randi 
variance a} and p is the correlation coefficient between (r) and 52(^)-Tl 
DOA’s were obtained by averaging over 100 Monte Carlo runs. The nun 
the number of virtual subarrays, the total array size and the particulars 
are described in the captions of figures and tables. The SNR indicated ir 
to the value at the input of the sensor element. For spectral MUSIC, the 
was conducted in steps of 0.002°. 

Figure 1 shows variation of the norm of the effective error vector ii 
steering vector (see (20)). Note that as the spacing between the eleme 
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Fable 1. Performance of the Root-MUSIC applied to transformed UCA data as a function of the cir 
:umferential spacing between the elements (without smoothing) (L = 30, A = 100, DOA's = 0° and 7° 
5NR = 3dB,p = 0.0). 

Circumferential _ MSE in DOA estimate (deg^) _ 

spacing between Spectral MUSIC applied Root-MUSIC applied to 

be elements to UCA data _ transformed UCA data _ 

[d) Simulation Evaluation of (14) 


132k 

0.2275 

0.07972 

0.08477 

134k 

0.1570 

0.06487 

0.05793 

136k 

0.1104 

0.06454 

0.05382 

13Sk 

0.07045 

0.04952 

0.03762 

3.40A. 

0.04468 

0.03284 

0.02803 

3.42A. 

0.03259 

0.03140 

0.02736 

).44A. 

0.02222 

0.04897 

0.03515 

D.46A. 

0.01627 

0.08968 

0.07188 

3.48A. 

0.01116 

0.30931 

0.25090 


the results. The performance of the spectral MUSIC improves as the spacing increases 
This is because of the increased aperture that we get with increasing value of d. Thi 
performance of the Root-MUSIC is better than that of the spectral MUSIC when th( 
spacing between the elements is less than 0.42A.. But, as the spacing is increased furthei 
the Root-MUSIC performance degrades and becomes worse as d approaches A/2. This i 
because, at larger spacings, the error in the transformed steering vector is quite large (se^ 
figure 1). The simulation results agree reasonably well with the theoretical values obtained 
from the numerical evaluation of (14) when d is less than 0.46A. For larger values of d 
however, the difference between the two increases since the theoretical expression (14) i 
less accurate when the term AR/7 becomes large. 

To see if the smoothing reduces the effect of the error introduced due to the approximatioi 
(cf. (3)) in the case of both uncorrelated and correlated scenarios, and to evaluate th 
utility of the theoretical result (18) and (19), we applied the Root-MUSIC with forwan 
smoothing to the covariance matrix (Rj-)w* The result so obtained from this is referred h 
as the asymptotic performance from the algorithm. Table 2 gives this result along with th 
theoretical values evaluated from (18) and (19) for various values of subarrays. Note fror 
the results of table 2 that the MSE is maximum for = 1 (no smoothing) and it drop 
significantly with smoothing. Consider the results shown in the table for uncorrealtei 
source scenario. As the smoothing is increased beyond K = S, the performance start 
deteriorating because of the reduction in the aperture. The results predicted from (19) ar 
not identical to those obtained from the algorithm since the theoretical expressions ar 
accurate for small values of errors and Aa(0), the error in the present case, is not small a 



UCA data tor two-source case, = ju, a = u.^oa, ukja s = u ana sj 


Number of 
virtual 
subarrays, 
{K) 


Asymptotic performance of the Root-M 
with forward smoothing (A0? in de^ 

P = 

= 0.0 


From the algorithm 

Evaluation of (19) 

From the algor 

1 

0.4852 X 10"' 

0.5359 X 10-' 

0.4852 X 10 

2 

0.1110 X 10-‘ 

0.1148 X 10-' 

0.2929 X 10 

3 

0.2002 X 10'* 

0.2018 X 10-' 

0.4686 X 10 

4 

0.4491 X 10-2 

0.4684 X 10-2 

0.2819 X 10 

5 

0.3120 X 10-2 

0.3594 X 10-2 

0.7611 X 10 

6 

0.5329 X 10-2 

0.5958 X 10-2 

0.9496 X 10 

7 

0.2543 X 10-2 

0.2624 X 10-2 

0.2612 X 10 

8 

0.1255 X 10-2 

0.1524 X 10-2 

0.5241 X 10 

9 

0.2663 X 10-2 

0.3085 X 10-2 

0.7199 X 10 

10 

0.3612 X 10-2 

0.4167 X 10-2 

0.8421 X 10 

11 

0.5708 X 10-2 

0.6524 X 10-2 

0.1081 X 10 

12 

0.1155 X 10-2 

0.1289 X 10-2 

0.1653 X 10 


the case of correlated sources are also given in table 2. Since the 
real, these results hold for both the forward and forward-backwan 
improves the performance in this case too. However, the improv 
with uncorrelated sources. 

To see the differential impact of forward and forward-backwar 
ence of highly correlated and closely spaced sources with finite d< 
nario with p = 0.95eJ^/^ and N (number of snapshots) = 100, ke 
nearly equal to the beamwidth of the UCA, and evaluated the h 
Table 3 gives the simulation results and the theoretical values prt 
Note that the forward-backward smoothing yields much superior p 


Table 3. Finite data performance of the Root-MUSIC with smoothing appl: 
with omni directional elements for correlated sources. (L = 50, d = 0.34X, 
10°, p = 0.95e^’^l\ SNR =3 dB) 


Number of 
virtual 
subarrays 
{K) 


MSE in DOA estimate (deg^) 

Spectral MUSIC 
applied to 

UCA data 

Root-MUSIC with smoothing applie 

Forward-Backward smoothing 


Simulation 

Evaluation of (23) 

Simu] 

1 

0.02485 

0.003313 

0.003058 

0.01 

2 


0.003419 

0.003223 

0.01 

3 


0.003568 

0.003449 

0.01 

4 


0.004373 

0.004053 

0.01 

5 


0.007673 

0.007213 

0.02 

6 


0.01883 

0.01789 

0.05 

7 


0.01569 

0.01582 

0.04 

8 


0.01462 

0.01556 

0.03 

9 


0.01817 

0.02131 

0.05 



;o the forward smoothing only. Further, the Root-MUSIC with forward-backward smooth- 
ng applied to the transformed UCA data performs better than the spectral MUSIC appliec 
;o the UCA data for all values of K (and much better at lower values of J^T). When K 
ncreases, the aperture comes down and the performance will start degrading when the 
iperture effect becomes predominant. The difference between the simulation results anc 
he predicted values (particularly in the forward smoothing case) can be attributed to the 
fact that the smoothing reduces the effect of finite data perturbations in addition to reducing 
he correlation among the sources, and hence, the actual MSE will be less than the value 
jiven by (14) (see Reddy & Reddy 1996a for discussion on this issue). 

5. Conclusions 

Fhis paper extends the spatial smoothing to UCA and analyzes the DOA estimation per 
formance of the Root-MUSIC with smoothing applied to the transformed UCA data. It h 
shown that the smoothing helps in reducing the effect of the errors that arise while extend 
ing the Root-MUSIC to UCA, in addition to reducing the correlation among the source: 
and the Effect of noise perturbations due to finite data. 
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Abstract. In this paper, we present a method for region-of-interest (ROI 
tomography using noisy projections. A wavelet decomposition down to th 
coarsest level is done on the noisy signal. The signal at various levels is estimate! 
using a Wiener filter. By assuming that the projections are 1// processes, th^ 
Wiener filtering reduces to a scalar multiplication. Using the Wiener filter an( 
regularity property of the wavelets, we combine the estimation and localisatioi 
of the noisy projections for ROI imaging. Experimental results are shown oi 
Shepp-Logan phantom and actual CT images. The validity of the 1 // mode 
for projections of real life images is also shown. 

Keywords. Region-of-interest reconstruction; noisy projection; fractal mod 
els; Wiener filtering. 


1. Introduction 

In computer-aided tomography, the objective is to reconstruct the cross-section of an objec 
from measurements that are strip integrals of some property of the object. In transmissioi 
X-ray tomography, the measurements consist of integrals (projections) of the attenuatioi 
coefficient /i(x, y) of the object along strips which represent the path of the X-rays througl 
the object. A popular technique for image reconstruction is the Filtered Back Projectio; 
(FBP) (Kak & Roberts 1986). 

In this work, we are concerned with the problem of reconstructing a portion of th 
image from noisy projections. The discussions in this paper are restricted to parallel bear 
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tomography. Reconstruction of only a portion of the cross-sectio 
of-interest tomography or local tomography) leads to reduction i 
to the patient and savings in computation (compared to FBP of tf 
However, the problem of local tomography is complicated b 
uniquely solvable in even dimensions (Natterer 1986). Most of i 
ducted in 2-D and the reconstruction formula becomes globally 
integrals of the object. Using wavelets is one method of solving t 
The FBP algorithm does not produce satisfactory reconstructioi 
The work reported here combines the ideas of multiresolution 
estimate of the projections using wavelets. We use a 1// model 1 
Wiener filter for the estimation. It is further shown that this algor 
case of ROI tomography. 

2. Preliminaries 

2.1 Radon transform and its inversion 

The two-dimensional Radon transform, * f (r, 6*), of a function / (x 
parameterized by (r, 9) such that r = x cos 0 -f y sin 0 is defined £ 

f(r, 0)= / /(x, y)5(r — X COS0 — y sin0)dxdy. 

The one-dimensional function fir, 9) is termed as the projectior 
9. The back-projection operator denoted as B is defined as (Kak < 

hsix, y) = Bh(r, ^) = f /i(xeos0-f y sin0, (r)d0, 
Jo 

Then it can be shown that 

f(x,y) = BJ^-\\co\mr,e)]. 

Equation (3) is the filtered back-projection implementation of ti 
form. It means that we filter each projection by a |a)| filter and thei 
projections. 

2.1a Nonlocality of the Radon inversion: The Radon inversio 
sented as follows (Olsen 1996): 

fix, y) = 7^-'f(r, 9) = BUr^ f(r, 9), 
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The Hilbert transform imposes a discontinuity on the Fourier transform of any functioi 
whose average value is not zero and also discontinuities on the higher derivatives whicl 
are not zero at the origin (Olsen 1996). This is because 

^'Hfit) = T (fit) * = T sign(a;)F(a>). (6 

Because of the discontinuities in the frequency domain, there is a spreading of the suppor 
of the function in the time domain (Olsen & DeStefano, 1994). However, this spreading 
of the function’s support will not occur if the function’s Fourier transform and the highe 
derivatives of it are zero at the origin (Olsen & DeStefano, 1994). If the Fourier transform o 
higher derivatives of it have zeroes at the origin, it implies that the function has a numbe 
of zero moments. Hence functions with arbitrarily high zero moments will have thei 
support unchanged by the Hilbert transform. It is in this context that the wavelet transfom 
is used. Wavelets are usually constructed with many zero moments (Daubechies 1988) 
Thus local properties of the high resolution components of a wavelet transform remaii 
local after applying the Hilbert transform. It has been shown in (Olsen & DeStefano, 1994 
Delany & Bresler, 1995) that using wavelets which have several zero moments, the wavele 
transformed filtered functions will have essential compact support. 

2.2 Wavelets 

2.2a Continuous wavelet transform of 1-D signals: The continuous wavelet transfom 
^fif)ia, b) of a signal fit) is defined as (Chui 1992): 

/ oo Z't — b\ 

^|ar‘/2/(0 jdr, (7 

where a € b e are the dilation and translation, respectively, of a single wavele 
function j/ it) called the mother wavelet, i{r denotes complex conjugation. A very importan 
property for wavelets is the regularity, which means the smoothness of the wavelet function 
Regularity R of the mother wavelet V^(r) is defined as (Aldroubi, 1996): 

/ +00 

firit)dt=0, (8 

-oo 

where r is an integer such that 0 < r < R. 

Thus a mother wavelet with regularity R has R — l vanishing moments or ^ ico) (Fourie 
transform of if it)) has a zero of order R at origin, because (Chui 1992): 


Given a function /, it is possible to obtain the low resolution si^ 
scale j in the following manner (Daubechies 1992): 

= (h{.) * 


and 

d^\k) = ig{.)*d^j+^\2k)), 

where * represents 1-D convolution, h and g are two filters assi 
finer scale approximation coefficients (/) can be synthesize 
scale approximation and detail coefficients, and as 

^ - k)f^j\n) + J^g{2n - k)d^ 

n n 

where ^ and | are defined as h(n) = hi-n) and g(n) = gi-n). 
2.4 1 // processes 

A 1// process is a nonstationary random process with a power S] 


where ct is a constant and y is known as the spectral exponent. 

It can be shown that (Womell & Oppenheim 1992) 

where D is known as the fractal dimension. If we choose gamm 
then 1 < D < 2. Hence fractal dimension of the signal is not an ii 
which lies in (1,2) and the signal is self-similar. This type of signj 
are known fractal signals. Here we concentrate on fractal sign 

2.5 Statistical properties of discrete wavelet transform coefficie 

From the regularity condition of the mother wavelet, 'F (cu), the F( 
has a zero of order R at the origin of the frequency plane. It ensure: 
with mother wavelet essentially removes the low frequency comp 
detail coefficients. These detail coefficients again are wide-sense si 
even when the signal is non-stationary (Masry 1993). This show 
concentrated aroimd zero frequency. In our analysis, we have a 
stationary. However, in case of any non-stationarity present in the 
will effectively remove it. Consequently, processing with only d 
be simnler than nrocessine the oririnal sienal itself. We have ad( 
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re 1. ACF of highest level detail coefficients with (a) Haar wavelet, (b) with 
echies’ wavelet. 


ot only makes detail signals WSS but almost decorrelates them also. Correlation 
3f detail coefficients can be obtained by calculating autocorrelation at a given 
over scale to scale. Denoting detail coefficients at scale m (resolution level 2^) 
^re m and n are the dilation and translation indicies respectively, it can be shown 
lan 1997) 

E[dj^dl] ~ 0{\ 2'"k - 2"/ (15) 

tion gives the correlation across scales, as well as at a particular scale. Conse- 
ith higher regularity mother wavelets it is possible to achieve fast decay of cor- 
)efficients. This property is illustrated in figure 2.5. Furthermore, it can be shown 
Jtail coefficients are mutually uncorrelated and their variance decreases geomet- 
h scale m and at each scale the variance, 0 * 2 ^, can be expressed as (Womell 1993) 

(16) 


er filtering 

assume that the original signal itself is a Gaussian process and corrupted by 
ndent additive white Gaussian noise (WGN). Clearly Wiener filtering gives an 
itimate of the uncorrupted signal since it minimises the mean square error in the 
case. Instead of doing Wiener filtering directly on the signal itself, we decompose 
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that wavelet decomposition whitens the 1// process. Complexit 
greatly reduced by adopting this technique (Womell & Oppenh( 
Let and x” be the detail coefficients of the corrupt s 

the original uncorrupt signal at a scale m (level 2'"). Hence, r™ 
doing Wiener filtering at each scale and the problem is to find 
response h'^ so that for the input the output will be an optimal 
by x'” the optimal estimate of x^ we have. 


-foo 


+00 




C*K= E E '"K-k. 

k— — OQ k=—OQ 


where ★ denotes linear convolution. 

The optimal filter impulse response at a given scale m is ] 


h 


m 

n 




S[nl 


where a and ^ are signal parameters, cr^ is the variance of the 
defined by the equation ct|,„ = where is the variant 

at resolution m. 

From the above equation we see that Wiener filter reduces i 
Here Wiener filtering is done only on detail coefficients becai 
of the energy is concentrated at low frequencies and conseque 
additive white Gaussian noise. 


4. ROI reconstruction 

We now consider the problem of ROI reconstruction when the in 
It is shown that projections at a different angles can be considei 
different values of y. 

Data Collection: The non-locality property of Radon transf 
reconstruction from projections only over that region. It requ 
ROI. However, projections can be taken for the entire object (wi 
at a coarse resolution) at a low number of angles and the val 
the missing angles interpolated suitably (Olsen & DeStefano 1 
jections at an increased number of angles are taken. A comprel 
data collection strategy can be found in (Delany & Bresler 1 
1994). 
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typical projections against the resolution on a log-log graph. T] 
in self-evident. 

The Wiener filtering process involves decomposition of the 
the length of the signal corresponds to the filter length. As ai 
Daubechies filter is used the signal length at the lowest level wil 
padding the signal with suitable number of zeros so that the sigi 
length X 2^ where A: is an integer, denoting the resolution level 
a signal of length 256, we zero-pad it upto 320. At level 5, the 
the filter length. The reconstructed signal, after an inverse wavel< 
estimate of the original signal. 


Estimating fractal signal parameters: Maximum likelihood (^ 
fractal signals was introduced by Womell & Oppenheim (1992 
decomposition technique is used as a tool and wavelet filter 
decorrelate the 1 // signal but provide a power law variance s 
also (figure 3). Here it is assumed that the original 1// signal 
additive white Gaussian noise and we estimate fractal dimension 
using Womell-Oppenheim ML estimation algorithm (Womell < 

4.1 Algorithm for ROI reconstruction 

Figure 4 outlines the entire filtering and reconstmction proces; 
should be noted that the diagram is simply for the purpose o; 
shown do not reflect the values we used in our simulation. N 
projections are sparsely sampled in the angular variable. In the 
N/2 full exposure projections and A/2 reduced exposure prc 
means that there are N projections over the ROI and A/2 away 
explained in the steps below. 

(1) A wavelet decomposition is done on each full exposure pro 
equal to filter length. The detail coefficients are estimated i 
adding only the detail coefficients corresponding to the prc 
obtain the estimate of the projections over the entire object, 

(2) A wavelet decomposition is done on the reduced exposure p: 
termediary angles and of smaller length. At the lowest possi 
is equal to the filter length. The low resolution coefficients c 
processed as they are not localised. Thus before doing a \ 
reduced exposure detail coefficients, we must obtain the 1 


Region-of-interest reconstruction from noisy projections 

'ull esqposure projections of length 40 
Leduced exposure projections of length 20 
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re 4. Simple schematic of the algorithm (the values tire simply for the purpose of 
nation; they do not reflect the values used in simulations), showing the wavelet decom- 
on, Wiener filtering and wavelet synthesis. Numerals denote signal lengths at various 
s. h and g are of length 10. 


using the detail coefficients available for the ROI. Note that the detail coefficients 
larticular resolution will be smaller in length than the low resolution coefficients, 
former is obtained only for the ROI. Thus the addition of the detail coefficients 
De only over the ROI. 

rojection data for all the angles is now obtained, with the detail information 
only to the ROI. 

)f the result obtained in step (4) results in an object reconstruction which has a 
esolution over the ROI and a low resolution elsewhere. 


)ping the synthesis of wavelet coefficients at a particular scale, we can obtain 
various resolutions. This can also be done sequentially, and the radiologist can 





















of the phantom image consists of 128 parallel rays while each projection of the 
images consist of 193 (approx. 128^2) rays. The ROI is defined to lie at the 
approximately covers l/4th of the image. The ROI can be shifted to some other 
hange of co-ordinates. We have used Daubechies’ wavelets of length 10 (/z is of 
(Daubechies 1992). 

5 5j, the image is reconstructed from the full projection data at 128 angles with- 
Dise. The effect of noise on the data when no estimation techniques are used 
in figure 5a, d and g, for 15, 10 and 5 dB. The result of using the filter with 
chies wavelet is shown in figure 5b, e and h. It is apparent that noise has been 
^ removed. The smoothing of the edges can be explained by looking at the corre- 
Drojection. A typical projection and the output after filtering is shown in figure 6. 
frequency details of the projections are removed as the filtering is unable to 
ite the high frequency signal from the noise. For purposes of comparison, we 
structed using the Haar basis. The presence of circular rings can be accounted 
ig the corresponding reconstructed projection. The estimated projection con- 
p edges. On back-projection, the sharp edges lead to rings in the image. The 
5S in the estimation using the Haar wavelet is because of several factors. The 
L process assumes that the wavelet coefficients are decorrelated from scale to 
are 1 shows the autocorrelation function with both the Haar and Daubechies’ 
The decorrelation in the second case is greater than in the first. Thus our as- 
is more correctly adhered to in the second case. As has been shown before, the 
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ure 8. Reconstruction of the images from noisy projections at 15 dB SNR. The columns 
b, c) show the original image, noisy image and noise removed image respectively. The 
^s correspond to images A, B, C respectively. 

'W good is the 1 // model? 

re algorithm described thus far hinges on the assumption that the projection of 
je at each angle is a fractal signal. To test our assumption, we considered a few 
edical images obtained from medical texts. Here we show the results on three of 
;es we collected (figure 8). We refer to the images as A, B and C. 

scan of the large frontal lobe. 

ital CT scan in the investigation of the cortical areas (frontal, temporal, parietal), 
ir, parasellar and third ventricle region. 

Dminal scan. 

9 shows the linear trend in the detail coefficients across scales. The next two 
ihow the noise removal and reconstruction usins this fractal model. The results 
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when projections are noisy For tWs 0 ^' ROI reconstn 

projections and filtered the noise ^ a 

output as it is done on the detailed Wiener filter. This filter 

position on the projections nese dei 
u^ed this property for ROI relon^^ 
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>nal invariance of two-level group codes over dihedral 
cyclic groups 
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Abstract, Phase-rotational invariance properties for two-level constructed, 
(using a binary code and a code over a residue class integer ring as component 
codes) Euclidean space codes (signal sets) in two and four dimensions are 
discussed. The label codes are group codes over dihedral and dicyclic groups 
respectively. A set of necessary and sufficient conditions on the component 
codes is obtained for the resulting signal sets to be rotationally invariant to 
several phase angles. 

Keywords. Multilevel codes; group codes; dihedral groups; coded 
modulation. 

duction 

mown (Viterbi & Omura 1979; Bendetto etal 1987) that digital communication 
Ltive White Gaussian Noise (AWGN) channel can be modelled as transmission 
from a finite set of points, called signal set, of a finite dimensional vector space 
laximum Likelihood soft decoding then becomes choosing the closest point in 
set, in the sense of Euclidean distance, from the received point in the space, 
ibility of error performance to a large extent is dominated by the minimum of 
ise distances of the signal points. The problem of signal set design for AWGN 
len is choosing a specified number of points in a space of specified dimensions 

minimum rliot^mr-A io mciYimnm nrvcciKlf> 
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where, dE(ci,tf) denotes the squared Euclidean distance between c 
identity element of G. If G and S have the same number of elemei 
of S can be labelled with the elements of G, and such a labelling sal 
(1) is referred as a matched labelling (Loeliger 1991, 1992). A sig 
group has the property that the Euclidean distance distribution of th 
set from any point is same, known as Uniform Error property (UEP) 
dimension N matched to a group G and /x is a matched labelling, the 
mapping 

/x": G" -> S" given by gi,.. •, gn-i) = if^(go), l^igi), ■ 

gives a signal set in Nn dimensions, called the signal space code and 
group C. Forney (1991) has shown that such signal space codes are sp 
class of codes known as geometrically uniform codes which have 1 
label code of the signal space code. 

Recently ‘multilevel constructions’ have been reported using binar 
codes and suitable mapping of coded bits onto a signal set of small dii 
a phase shift keying (PSK) signal set (Cusack 1984; Sayegh 1986; I 
Kschischang et al 1989; Kasami et al 1991; Calderbank & Sesha 
Caire 1994; Garello & Bendetto 1995; Imai & Hirakawa 1997). M 
has attracted wide spread research because of their amenability foi 
multistage decoding (Calderbank 1989; Takata era/1993; Kofman et i 
block code of L levels uses L block codes, called component codes, 
over finite alphabets of possibly different sizes. A signal set 5, calk 
dimension N, has O/Ci tni points, where m,-, / = 1,2,..., L, are the 
with each point labelled by an ordered L-tuple with one entry froir 
this labelling, a set of L codewords, one from each code, correspo 
dimensions, with each coordinate of L codewords choosing a point i: 
all such points corresponding to all possible combinations of codewc 
codes is the multilevel constructed signal space code or Euclidean s; 

This paper reports rotational invariance properties of coded sip 
two-level (L = 2) construction for a class of two-dimensional {N - 
(SPSK) and four-dimensional (iV = 4) signal sets. Four-dimensional 
studied by several authors (Welti & Lee 1974; Zetterberg & Branc 
Birdsall 1989; Visintin et al 1992). Gersho & Lawrence (1984) des( 
and implementation for a particular 2-bits per dimension four-dimens 
readily lends itself to simple encoding and decoding. For this encodinj 
gain in noise margin over conventional 16-point(two-dimensional) 
modulation (QAM) signalling. In our two-level construction the com] 
a binary code and a linear code over a residue class integer ring. The 


’he class of dihedral groups is defined by 

= e, r‘s = sr~‘, 0 < i < M, j = 0, 1} where e is the identity 

»f Dm- The group operation can be expressed as 
ind the inverse of an element is given by 

fhis group has 2M elements, where M is an arbitrary integer and r and 5 are called the 
generators of Dm- 

DEFINITION I 

Let S denote a unit circle in 2-dimensions. A matched labelling ix : Dm S, for a 2M- 
isymmetric PSK signal set matched to Dm is said to be an m-labelling, 0 < m < M - 1 
with angle of asymmetry (/>, —n/2M < <p < 7i/2M, and denoted by mL(i,, if 

i = 0,1,...,M-1, y =0,1,(/,M) = 1 

Definition 1 is general and includes Asymmetric PSK(APSK) signal sets. It is showi 
(Bali & Rajan 1997b) that for a given group code APSK performs better than SPSK unde 
certain conditions. In this paper, however, we will be considering only SPSK signal set 
i.e., 0 = 0 throughout. Observe that various values of the parameter m gives differen 
labellings of the SPSK signal set with m = 0 corresponding to labelling the signal point 
in natural order in the anticlockwise direction. Our results in this paper hold even if m C 
The class of dicyclic groups is defined by 

DC2M = 0 < i < 2M, j = 0, 1}. 

This group has AM elements, M being any positive integer and is generated by r and 
where, = e, the identity element. The group operation can be expressed as 

)(^!2^)2) — ^(ii+<2+)l(^y2—2i2))modulo2M^(yi-l-j2)w'od'ilo2 

and the inverse of an element is given by 

DEFINITION 2 (Bali & Rajan 1997b) 

Let S denote the unit sphere in 4 dimensions. The matched labelling we consider is tl 
subset of S (shown in figure 1 - the - 2:2 plane constitutes the first two dimensions ar 

_ IT., mmi^initicr tWoV 
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consisting of AM points and the matched labelling jx : DC 2 M is 

= (1 — 1) cos (Ictc/M) , (1 — /) sin {kn/M), 

I cos {kn/M ), — 1 sin (kn/M). 

It is routine calculation to check that the mappings of definitions 1 and : 
the conditions for matched labelling given in (1). 


3. Two-level group codes and their characterization 

The block diagram of a two-level block-coded modulation is shown in figu 
are length n codes over alphabets X = {xi, a: 2 , X 3 , X 4 }, (m 1 =4) and Y = 
2). Figure 2b shows a labelling of S consisting of eight points on the circh 
of X and Y. For codewords a = (ao ,..., a„_i) € Ci and b = (bo, ■ ■ 
each pair (at, bi),i =0, 1,..., n — 1, selects a point in 5, and the pair ( 
point in 2n dimensions. The collection of all such points in 2n dimension 
to all possible pairs of codewords constitute the two-level block coded n 
(signal set) or signal space code. This paper deals with X and Y being Z 2 
class integers modulo 2 and modulo X respectively, where X = M for cod 
X = 2M for codes over DC2m- 

DEFINITION 3 

Let Cs and Cr be length n codes over respectively, Z 2 = {0,1} and Zx = { 
and a = (ao,..., an-i) ^ Cs,b = (bo,..., b„-i) e Cr- Let Ca,b denote 


Rotational invariance of two-level group codes A{ 



Figure 2. (a) Block diagram of a two-level block coded modulation, (b) Labelling of an 

8-PSK signal set with X and Y. 

Bigleiri & Caire (1994) have studied the Geometrical Uniformity properties of signa 
space codes obtained by L-level construction. The component codes used are L linea 
binary codes with the 2^-ary SPSK as the basic signal set. The points of the PSK signa 
set can be designated with either the cyclic group with 2^ elements or with the dihedra 
group with 2^ elements. They derive conditions under which the resulting multilevel cod 
is a group code over the cyclic or dihedral group. 

In the rest of the section we give the necessary and sufficient conditions on the componen 
codes of the two-level construction shown in figure 2a to result in a group code over Da 
and DC 2 M- 

Theorem 1 (Bali & Rajan 1997a, 1998). The two-level code C = r^'s^^ is a group cod 
over Dm if and only if 
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obtained as 

Cr = Cl + 2C2+• ■ • • + 2 ^ 

then theorem 1 coincides with theorem 1 of Bigleiri & Caire (1994). 

Proof of theorem 1 and the Euclidean distance properties of the resulti 
asymmetric PSK signal space codes can be seen in Bali & Rajan (1998 

Example 1. If Cr = Z^, i.e., if Cr is the trivial code consisting of all 
over Zm, then for any binary linear code Cs, conditions of theorem 1 
hence the resulting code is a group code over Dm- 

Example 2. If Cs is a repetition code consisting of all zero and all om 
Cr is any linear code over Zm, then conditions of theorem 1 are satis! 
code is a group over Dm- 


Example 3. LetC^ = {000, 111} and = (000, 220,132, 312, 200,0: 
Z 4 .Then the resulting two-level code is a group code over D 4 and has 
codewords 


{(r°, r°, r°), (r^, r^, r°), (r^r^,r^), (r^,r\r^), 
(^2, ^0^ ^0^^ ^^0^ ^2^ ^0^^ ^^3^ ^3^ ^2^^ ^^1^ r^ r^), 

(5, 5, s), (r^s, r^s, s), (r^s, r^s, r^s), (r^s, r^s, r^s), 
(r^, s, s), (^, r^s, 5), (r^s, r^s, r^s), (r^s, r*s, r^ 5 )}. 


Observe that in this example, M is a power of 2 , but Cr is a non-dec< 


Example 4. Let Cr = (0000,1233,3211,0222,1011, 2022,3033,2: 
Cs = (0000, nil). For these codes Cs O 2Cr = (0000, 2022} C ( 
two-level code jg ^ group code over D 4 and the codewords are; 

( 1111 ), (l,r^,r^,r^), (r^,r^,r,r), (r, r^, r^, r^), 

(r, 1 , r, r), (r^,r^, 1 , 1 ), (r^, l,r^.r^), (r^, l,r^, r^), 

(s, s, s, s), (s, r^s, r^s, r^s), (r^s, r^s, rs, rs), 

(rs, r^s, r^s, Ps), {rs, s, rs, rs), {r^s, r^s, s, s), 

{r^s, s, r^s.r^s), {r^, 1 , r^, r^), 


The following theorem gives the necessary and sufficient conditions 
codes of the two-level construction shown in figure 2 a to result in a group 
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The proof is similar to that of theorem 1 and is omitted. We list below a few classt 
group codes over DCx- 

Example 5. If Cr = i.e., if Cr is the trivial code consisting of all possible «-tu 

over Z 2 X, then for any binary linear code Cj, conditions of theorem 1 are satisfied 
hence the resulting code is a group code over DCx- 

Example 6. If Cs is a repetition code consisting of all zero and all one vectors only 
conditions of theorem 1 are satisfied for any linear code over Z 2 x containing all 
codeword, as Cr code. 

Example 7. Let Q be any binary code with the property Cs Q Cs Q Cs and Cr is 
parity check code over Z 2 x- Then the resulting two-level code is a group code over £ 

Notice that when all the three conditions of theorem I are satisfied for a set of compo 
codes then these can be used to construct both group codes over the dihedral group as 
as group codes over the dicyclic group. The resulting coding gain when used as dihe 
group code is given by 

Gd = dD/^sin^in/l^), 

and when used as dicyclic group code is given by 

where do and d^c denote the MSED of the signal space codes of dihedral and dic} 
group codes respectively and R denotes the rate of the code in bits per channel 
(Kschischang et al 1989). 

Theorems 3 and 4 below identify two situations for component codes satisfying 
ditions of both theorems 1 and 2; labelling with dicyclic groups gives better perform; 
compared to labelling with dihedral groups. 

Theorem 3. IfCs isa{n,k,ds) binary linear code with CsOCs c Q and Cr = Z 2 M 
then, the two-level code construction will give group code over both Dm <ind DC 2 M 
such a pair of component codes 

(i) If 2cos(nr/4M) > TiR), 

where r(i?) = sin(7r/2^/^)/sin(7r/2^) then, Gx > Gp. 

(ii) If R > 2.68 then, Gp > Gx- 


52 


J Bali and B S Rajan 


Theorem 4. IfCs is a (n, n/2, 2) binary code (n even) with gef 

ril0000...00“| 

0 0 1 1 0 0 ... 0 0 


Loooooo... iij 

and Cr is {n,n ~ 2) 2M-ary parity check code then, Gx - 

sin(7r/2^/^)/ sin(7r/2^). 

Proof. It is easy to verify that the component codes in the stateir 
conditions of both theorems 1 and 2. The rest of the proof is si 
theorem 3. 

Example 8. Consider the case n = S, M = 2 in theorem 4. 
bits/symbol, 

Gd = 2sin2(;r/8)/sin2(;r/22-25) = 0.7782735 (-1.0 
2sin2(7r/4)/sin2(7r/22-25/2) = 1.0171907 (0.0740239 dB) > ( 

4. Phase rotational invariance 

For L-level codes with binary component codes phase invaria 
reported (Kasami et al 1994a). In this section we discuss phase i 
erties of two-level group codes over Dm and DC2m- 
For mL of a 2M-SPSK signal set, a codeword 

(^yos^o^ ^ ..., 6 C, 

is mapped onto the point 

^^p{V^[xo((2m+l)7r/M)+yo27t/M]}^ . . ., 

in 2n-dimensional space. There is a one-one correspondence be 
dimensional points given by (3) and (4). The code is said to be i 
angle if whenever (4) is a signal point for a codeword in C, thei 
to 


^^„lV^lxo((2m+l)^)+yo27r/M+0]} 



'litjM), rotations where k divides M, iff the all k vector (k, k, ..., k) e Cr, 
' 27 tlM)y rotations where k and M are relatively prime, iff the all-1 vector (I, 1, ..., 
€ Cry 

IM rotations iff all-1 vector is present in Cr and Cs> 


>iM/k = X. Then replacing 3 ^/ by 3 ;/ + X, / = 0, 1,..., /i — 1, in (3) corresponds 
6 in ( 4 ) getting replaced by kljtfM, and the converse. 

k and M are relatively prime, then klrc/M rotations can be obtained by k successive 
r/M rotations and In jM rotations can be obtained by u successive kljtjM rotations 
here u is given by uk + vM = 1 (Bezout’s Theorem). Hence it is sufficient to 
insider In/M rotations only for which all-1 vector (1, 1,..., 1 ) should be in Q, 
hich follows from (i) with X = M. 

ippose all-1 vectors are present in both Cr and Q. Presence of all-1 vector in Q 
larantees rotational invariance by ( 2 m-|-l) 7 r/M = m27r/M-f7r/M.Thepresenceof 
I -1 vector in Cr guarantees rotational invariance by all multiples of In jM, including 
mln/M, Clearly, rotational invariance for both mlTtjM + jr/M and —mlix/M 
iplies rotational invariance for n/M. The converse is straightforward. □ 


a 4 n-dimensional signal set obtained using a two level group code over DC 2 m> ^ 
3rd 


ped onto the point 

(1 - a:o) xo ... 

(1 - Xi) xi ... 

(1 - x„_i) xn-x 


(5) 


( 6 ) 


limensional space. There is a one-one correspondence between codewords of the 
code given by (5) and 4n-dimensional points given by ( 6 ). The code is said to be 
nally invariant to an angle 0, if whenever ( 6 ) is a signal point for a codeword in C, 
le codeword corresponding to 

(1 - xo) exp''^[>'o^/^+^^\ xo ... 

(1 - Xi) exptV=T[yi7r/M+e]i^ ^ _ _ 

(1 - x„_i) x„_i 

leword in C. The following theorem give the conditions on the component codes for 
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Theorem 6. A two-level group code, C DCim invariant ti 

(i) k{n/M), rotations where k divides M, iff the all k vector {k,k,,.. ,k) ' 

(ii) k{ 7 t/M), rotations where kand Mare relatively prime, iff the all-1 vector 

Cr, 


Notice that the minimum angle of rotational invariance in this case is 7 r/i 
2M-SPSK it is In/M, Moreover, the G code plays no role in the condition 
invariance, i.e., the minimal angle of rotational invariance is completely dete 
Cr code. 

5. Conclusion 

Rotational invariance properties of signal space codes obtained by group < 
hedral groups and dicyclic groups are discussed. The group codes are th( 
by two-level construction from a binary code and a code over residue clas 
The signal space codes obtained from group codes over dihedral group C( 
dimensions and those obtained from dicyclic group codes are in An dimensi 
the length of the code. Theorem 3 includes conditions under which the code 
to the minimum angle for dihedral group codes. If the component codes s; 
ditions of the theorem then the effect of phase ambiguity in carrier recovei 
of InIM can be nullified by two-stage differential encoding as follows: 
coder (and before the modulator) the sequence of symbols to be transmitter 
(..., yiXi, yi+ixi^u • • •) where (..., x/, Xf+i,...) is the 

binary sequence and , yt-i, y/, y/+i,...) is the corresponding sequenc 
First differentially encode the binary sequence {x/} resulting in the seque 
After this the original sequence becomes (..., yz-ix^^_p y/xf, y/-i-ix^_^_j,.. 
above sequence can be seen as the interleaved version of two subsequences 
subsequence, say Sq = {yyxj}, consists of those symbols with xj = 0 and t 
Si = {yjXj], consists of those symbols with xj = 1. Now differentially e 
M, the [yj] parts of the each subsequence. 

It is straightforward to verify that by two stage differential decoding pn 
reverse order of encoding, the effects of phase rotations by multiple of th 
can be removed. 

Theorem 4 includes conditions under which the codes are invariant to the r 
for dicyclic group codes. If the component codes satisfy the conditions of th 
the effect of phase ambiguity in carrier recovery in multiples of n/M can 
differential encoding as in the dihedral group code case without the first s 
i.e., differential encoding of binary symbols is not needed for due to rotatioi 
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ipplication to Wigner-Ville distribution: A group delay 
ipproach 

S V NARASIMHAN *, EIPLOTKIN and M N S SWAMY 

Centre for Communication and Signal Processing, Department of Electrical and 
Computer Engineering, Concordia University, 1455, De Maisonneuve Blvd 
West, Montreal, H3G IMS, Canada 

^Present address: Aerospace Electronics Division, National Aerospace Labs, 
Bangalore, 560017, India 
e-mail: svnarsim@nalsic.emet.in 

Abstract. In this paper, a method of estimating the power spectrum of t 
complex signal based on the Group Delay function (GD) is proposed and also 
applied to the Wigner-Ville Distribution (WVD) to reduce the ripple effec 
due to the truncation of the autocorrelation sequence. The proposed methoc 
is realised by the GD for a complex signal and the modified GD concept 
This extends the performance advantages of the modified GD applicable to £ 
real signal, to a complex one. Further, its application to WVD, reduces th< 
truncation/ripple effect without sacrificing any frequency resolution, as nc 
common window function is used. Significant improvement in performance 
in terms of reduction in variance without any compromise on resolution anc 
higher noise immunity, has been found over those of the periodogram anc 
windowed WVD. 

Keywords. Complex signals; power spectrum estimation; group delay ap 
proach; Wigner-Ville distribution; Gibb’s ripple effect. 


1. Introduction 

aenerally, spectral estimation aims at extracting information about the system, from it 
Dbserved output (in the absence of the input), and a signal, when it is associated witl 
noise. A good spectral estimator, for a given length of data, provides an estimate whicl 
[las as high a resolution and as low a bias and variance as possible. These factors depem 
Dn the signal scenario (data length, window and the associated noise). The variance o 
the basic periodogram is large, though it has good resolution, low bias and good signs 


lobe of the spectrum of the window and this results in a poor resoli 
the cost of resolution, the averaged PSE provides a lower variance, c 
number of segments averaged and the window used, for a given length 
based parametric methods (Kay 1988), even for a relatively short d 
high resolution and low variance. But this is valid only when the sign; 
assumed model and the signal to noise ratio is high. 

Though the power spectral density (PSD) based on the Fourier tr; 
spectral information, it does not explicitly provide a time reckoning 
this, the spectral content is the same for ever. To analyse nonstatioi 
time FT, which gives spectral information about the data within a win 
was introduced and this involves a tradeoff between the time localiza 
resolution. To alleviate this, the Wigner-Ville distribution (WVD), 
representation (TFR) was introduced (Cohen 1989). The WVD, at an> 
FT of the instantaneous autocorrelation (lACR) sequence of infinite laj 
theoretically, it has infinite resolution both in time and frequency. H( 
it is the Pseudo WVD (PWVD) that is computed which considers IA< 
number of lags. In PWVD, in order to overcome the abrupt trunca) 
ripples (Gibb’s effect) in the TFR along the frequency axis, the lA 
a common window function and for a given lag length, this deterioi 
resolution. 

The group delay function (GD), the negative derivative of phase 
provides improved resolution over that of the PSE and facilitates 
(Yegnanarayana et al 1984). However, spectral estimation based on 
PSE, also suffers from smearing and the reduction in variance is onl 
resolution. Recently, a modification for the GD, which not only presei 
of the PSE, in particular, the good frequency resolution, but also reduce; 
estimate, has been proposed (Murthy & Yegnanarayana 1991; Yegnar 
1992). This basically removes the zeros close to the unit circle, without i 
and zeros of the system or the signal (assuming the zeros if any are c 
White noise or spectral ripples introduce zeros close to the unit circle, 
the output of a system driven by a white noise or a signal with its assoc: 
signal spectrum associated with ripples due to signal truncation, the va 
zeros close to the unit circle, irrespective of their origin. Consequent 
effect of only these zeros results in reduced variance without compron 
resolution. 

In this paper, modification for the GD has been applied to a con 
achieved by considering the GD for a complex signal (Reddy & Rao ] 
WVD, the lACR being an analytic complex signal, the modified grouf 
a complex signal (GDCM) approach has been tailored to remove the ri 
applying any common window, to achieve preservation of the freque 
better noise immunity. 


l.\ Group delay function for a complex signal 

[f x(n) is a minimum phase signal, then (Reddy & Rao 1987) 

oo 

In |X(£t))l = ^[c/?(n) cos cure + c/(re) sin cure], (!’ 

n=0 

00 

0 (cu) = —CR{n) ^mcon + cj{n) coscon], (2' 

n—0 

where 9(co) = 0u(cu) + Ire 1(a)), 9y(a)) is the principal value of the phase, 9(a)) is the 
unwrapped phase and c(n) = CR(n) + jci (re) are the complex cepstral coefficients. R anc 
I refer to the real and imaginary parts. For a minimum phase signal, the logmagnitude 
spectrum and the phase are related by a single set of complex cepstral coefficients. The 
GD Tm (co) is given by 

Xm(o)) = —d9(a))lda) 

OO 

= "^[ncR(n) coscon + ncj(n) sin con]. (3 

n=Q 

CR(n) and cj (n) are derived from the magnitude and hence Xmio)) is called the magnitude 
GD for a complex signal (MGDC). Also, (3) is 

= (l/2)FT[(nc(n) - nc^(-n))l (4 

where * indicates conjugate operation. If nc(n) is conjugate symmetric, then 

tmico) = FT[nc(n)l (5 

For a mixed phase signal, the logmagnitude and phase are not related by a single se 
of cepstral coefficients and two different cepstral sequences ci(n) = ciR(n) + jc\i(n 
and C2(n) = C2R(n) + jc2i(n), for magnitude and phase respectively, are defined. c\(n 
and C 2 (n) are the complex cepstral coefficients of minimum phase signals derived fron 
magnitude and phase respectively. For this case, the magnitude group delay (MGDC) i 
the same as given by (3) (using ci(n) instead of c(n)); however, the Phase group dela; 
(PGDC) is different and is given by 
00 

Tpico) = "^[nc2R(n) COS con+ nc2i(n) sin con], (6 

To compute, rp(co), from ( 6 ), C 2 {n) has to be derived from the phase function 9{co 
which is an unwrapped one. But tpico) can be directly computed from the FTs X(co) nn 
F(a;) of the signals x(n) and yin) = nx(n), respectively, by 

rp(a)) = -[XR(a))YR(a)) + Xi(a))Yj(a))]/\X(a))\^ 


0 


PGDC can be used for computing the MGDC. 

2.2 The modified group delay for a complex signal 


The variance of a spectral estimate is due to the undesired fine structun 
capture the spectral envelope and discard the fine structure without ai 
The zeros close to the unit circle manifest as spikes in the GD. Furtl 
the spikes depend upon whether the zero is inside or outside the uni 
spikes which contribute significantly to the fine structure of the spectr 
cannot be removed by normal smoothing without the loss of resolutior 
suggested by Murthy & Yegnanarayana (Murthy & Yegnanarayana 19' 
& Murthy 1992) removes these spikes effectively. Presently, it is pro] 
modification to the MGD for complex signals (MGDC). 

Assuming that the signal under consideration x(n) is generated by 
system, driven by a complex white noise or has complex sinusoids v 
noise and, further, that its spectrum X(w) can be put in a rational fon 
denominator corresponds to the system or sinusoids and the numerator 
excitation or the associated noise, respectively. For this case, the MGE 

Xmico) = 

where and ZmD. (<^), are the MGDCs for N(a)) and D(co) 

is given by (7), where x{yi) ~ Xjfi (n) and y{n) = ymin) = nx 
the minimum phase equivalent of x{n). It can be shown that 
_ aN(.o^) ctpjM) 

%(&>)- |^(^)p |0(^)|2- 

aAf (to) and a© (o)) are the numerators of (7) for n (a>) and Xmoico) re 

will have large amplitude spikes due to the very small values of 
which are close to the unit circle, while this is not so with (m), as th 
well within the unit circle. Hence, in % (co), the effect of excitation or t 
(both have zeros near the unit circle) masks the system or the signal o 
assumed to be an all-pole one. The effect of these zeros could be redu 
(9), by |A/^(m)|^. Further, as the envelope of |7/(co)p is nearly 
features of imP (co) continue to exist, with \N{cl))\^ fluctuations superim 
the modified MGDC (MGDCM) tmoico), is 

Tmo(co) = Xm(co)\TSI ((JL>)\^. 

For the computation of MGDCM, an estimate of 771 (uj) p, |iV (ru) p is i 
to the MGDC. This is achieved by computing the ratio of the sign^ 
\X{o)) p, to the smoothed power spectrum of the signal, |X((u) p, obtain 
cepstral coefficient sequence. The values of \N(o))\^ around the zeros h 
so that it cancels small values in the denominator of the first term of 
estimate | iV (tt)) p should retain all the sharp fluctuations of the logmagn 


of GD from the magnitude spectram is different from that for a real signal. Further, fc 
a minimum phase signal, as the MGDC and PGDC are same, the above expression fc 
PGDC can be used for computing the MGDC. 

2.2 The modified group delay for a complex signal 


The variance of a spectral estimate is due to the undesired fine structure. It is of interest t 
capture the spectral envelope and discard the fine structure without affecting the forme: 
The zeros close to the unit circle manifest as spikes in the GD. Further, the polarity o 
the spikes depend upon whether the zero is inside or outside the unit circle. It is thes 
spikes which contribute significantly to the fine structure of the spectrum and their effec 
cannot be removed by normal smoothing without the loss of resolution. The modificatio 
suggested by Murthy & Yegnanarayana (Murthy & Yegnanarayana 1991; Yegnanarayan 
& Murthy 1992) removes these spikes effectively. Presently, it is proposed to apply thi 
modification to the MGD for complex signals (MGDC). 

Assuming that the signal under consideration x{n) is generated by a complex all-pol 
system, driven by a complex white noise or has complex sinusoids with complex whit 
noise and, further, that its spectrum X(w) can be put in a rational form N(co)/D(co), th 
denominator corresponds to the system or sinusoids and the numerator corresponds to th 
excitation or the associated noise, respectively. For this case, the MGDC is given by 

tmico) = r,„N(co) - Tmoico), (8 


where and Xmo, i^), are the MGDCs for Nico) and D(co) respectively. Alsc 

Xmico) is given by (7), where x{n) = Xm(n) and y(n) = y,n{n) = nxm(n), x,„(n) bein 
the minimum phase equivalent of x(n). It can be shown that 


Xmico) = 


ap^ico) aoico) 


(9 


liV(cu)|2 |D(cu)|2- 
aN(o)) andaoCtv) are the numerators of (7) for rmNixo) and Xmoixo) respectively. 
will have large amplitude spikes due to the very small values of |fV((y)|2, near the zero 
which are close to the unit circle, while this is not so with Xmo (tv), as the roots of D {co) ar 
well within the unit circle. Hence, in x^ {co), the effect of excitation or the associated nois 
(both have zeros near the unit circle) masks the system or the signal component which i 
assumed to be an all-pole one. The effect of these zeros could be reduced by multiplyin, 
Xmico), (9), by iA((y)|2. Further, as the envelope of |iV(ct;)|2 is nearly flat, the significar 
features of x^d (co) continue to exist, with | A/'((y) |2 fluctuations superimposed on it. Hence 
the modified MGDC (MGDCM) x,„oico), is 

^mo {co) = x„{co)\N{co)\^. 


(10 


For the computation of MGDCM, an estimate of A | (cu) |2, | iV(cu) |2 is required in additio: 
to the MGDC. This is achieved by computing the ratio of the signal power spectrun 
IX (ru) 12, to the smoothed power spectrum of the signal, | Y (<z>) |2, obtained by the truncatei 
cepstral coefficient sequence. The values of | A(a)) (2 around the zeros have to be preserve 


complex signal 


"or a signal x(r), WVD is defined as 


WAt,co) 


-j: 


x(/+ r/2 )x*(f — T/2)e 


(11 


where r(r) = [x(r + rj2)x*it — r/2)] is the instantaneous autocorrelation function ant 
1 = indicates conjugate operation. For computational purposes, it is necessary to weigh th< 
signal by a window before evaluating the WVD and this window slides along the time axis 
"or a window function, h(t), h{t) = 0 for |t| > T/2, the WVD of the windowed signal is 


Wxhit, CO) 





(12 


where Wh (t, co) is the WVD of the window function. This WVD of the windowed signal i: 
called Pseudo Wigner-Ville Distribution {PWVD), PW^(r, a>). The effect of the window 
s to smear the WVD along the frequency axis. For a real symmetrical window, the PWVE 
s the FT of the windowed function [x{t + T/2)x*(t — t/2)], the window being /i^(t/2) 
rhe window eats away the coirelation function at higher lags which results in poor spectra 
•esolution. 

The quadratic operation on the signal, causes the WVD to be bilinear and this intro 
luces crossterms for multicomponent signals (Flandrin 1984; Velez & Absher 1990). Th( 
irossterm appears midway between every two components of the signal. Its amplitude i: 
Droportional to the product of the two components’ amplitudes and it oscillates in time a 
i frequency equal to the frequency separation between them. The crossterm effect, whicl 
nakes the interpretation of the WVD difficult, can be reduced by smoothing the WVI 
dong the time axis. The smoothing process, in time for crossterms and in frequency for th( 
ag window, can be considered a two-dimensional convolution of the WVD of the signa 
vith that of the smoothing kernel, d>(t, n>) (Flandrin 1984). Using different smoothing 
cemels, a class of distributions known as Cohen’s class (Cohen 1989) can be realized 
rhe kernel determines the properties of the distribution (Jeong Sc Williams 1992). Tb 
properties, viz. marginality in frequency, instantaneous frequency and frequency support 
ire not satisfied for common windows due to smearing. 

The very definition of the WVD demands the signal to be sampled at twice the Nyquis 
sampling rate, otherwise it introduces aliasing in the frequency domain as the periodicit; 
is n instead of lit. This is overcome by using an analytic signal (Picone 1988) whicl 
necessitates further processing techniques to handle complex signals. 

In the computation of PWVD, for a given lag length, the lag windowing used to tak 
:are of the truncation effects eats away the higher lACR lags and reduces the frequenc; 
resolution relative to that of no windowing. Further, the very definition of the WVD, whici 
results in an aliasing in the frequency domain, demands the conversion of a real signt 
to a complex analytic signal. With the MGDCM, as the zeros close to the unit circle ar 
removed without applying any window, the loss in the frequency resolution which woul 
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the zeros due to the associated white noise, the proposed WVD indirec 
noise immunity. 

For the MGDC computation, the starting point is the signal undei 
mentioned in § 2.1. However, to apply the MGDC to the present TFR, tl 
along the time axis (SIACR) is the beginning point and the magnitude 
(1), at a particular time instant is obtained from the FT of the SIACR 
WVD slice estimate. This estimate, at a frequency, is supposed to be a po 
represents the power spectral density (PSD). However, due to the inevita 
rectangular window, the values may become negative. This is due to the 
convolved with the FT of the rectangular window. Presently, as the c 
MGDC T^Cm), involves logarithmic operation, ((1) and (3)), it is neces 
the FT of the SIACR is positive. This is achieved by raising the floor le 
scaling up the SIACR at zeroth lag. Further, in computing (<w) the eqi 
spectrum is computed from the PSD and the linearly weighted cepstral C( 
is made conjugate symmetric. 

According to (10), the MGDCM, _rmo(<^), is computed by using an ei 
and this estimate is [X (o))/X{(a)]^,X (o)) being the cepstrally smoothed 
by the truncated cepstral sequence of x(n). That is, 

iV(£u) = X{a>)/X{a) = 1 + [A(m)/X(m)]. 

Here, A(a)) represents the fluctuating part of X(co). By multiplying 
spectmm which is free from contribution due to input or the associatet 
improved resolution is obtained and this is different from X (&>). For a s 
spectral characteristic, in the GD tmico), the contribution is only due tc 
to get a Tmoica), which is free from fluctuations, Xm(co) has to be multi] 

= rm((u)|A(<u)|^. 

Presently, it is required to remove the ripple on the floor, which is ( 
spectral characteristic and the pedestal height is immaterial as it is not n 
The PSD which is free from the ringing/ripple effect and has better res 
from Xmoicx)), using (14) and by retracing the MGDC computation proce 
order ((3) and (1)). It is important to note that in obtaining the processet 
coefficient sequence has to be made conjugate symmetric. If only tl 
importance, then only the positive values of the Xmoito) above a certai 
considered. To achieve the improved WVD, these operations have to be 
sample. It is important to note that since no window is applied, the propc 
frequency domain properties better than the PWVD. 


sinusoids in complex white Gaussian noise. The AR process and the complex 
; are the same as those considered by (Yegnanarayana & Murthy 1992), except 
brmer is a complex version. 

R system considered has roots at (0.6321 + 70 . 7593 ) and (0.7569 + 70.6204). 
the same AR process which is associated with complex white Gaussian noise 
signal-to-noise ratio (SNR) equal to 0 dB is considered. For sinusoids plus white 


x(n) = + V20e^2^(0’5)« + „(„), (15) 

ered the signal. u(n) is the additive complex white Gaussian noise, and its variance 
id to get SNRs equal to 20 dB and 0 dB. For the AR process and for the sinusoids 
le, the data lengths used are 256 and 100 points respectively. For both the cases, 
iting the |lV(a))p estimate, the smoothed spectrum is obtained by considering the 
:epstral coefficients. 

jiformance of the MGDCM for the AR signal and for the sinusoids considered, 
= 0 dB, are shown in figures 1 and 2 respectively. Further, these results are com- 
;th those of the PSE. In figures 1 and 2, (a) and (c) correspond to individual 
5 ; (b) and (d) show the mean and variance of the estimate obtained by 50 reali- 
rhe individual estimate obtained by MGDCM is free from the fine structure due 
put and the associated noise. The variance plots shown also support this as the 
ance is very high compared to those of MGDCM. Further, the reduction in vari- 
ieved by the proposed method over the PSE is quantified by computing the Sum 
Variance Ratio (SSVR) given by SSVR = EVz(co)/i:Vp(w), where Vr{o)) 
»r the variance computed by the GD approach and Vp{co) for that by PSE. For 
process and for sinusoids, the SSVR has values of 0.1783 and 0.0438 respec- 
he mean plots are similar, except that the mean obtained by the PSE has still 
e fluctuations, even after averaging 50 estimates. For the clean AR process and 
inusoids with noise (SNR = 20 dB), the SSVR values are 0.7471 and 0.0126 re- 
y. The reduction in variance achieved is more significant for sinusoids than for 
irocess. 

moothing parameter, the number of cepstral coefficients used, determines the 
of the proposed estimators and to get a good variance reduction, this number 
e kept as small as possible. Compared to the PSE, the proposed estimator has a 
ntly low variance. Further, depending upon the signal and its associated noise, its 
ranges only from 70% to 2% of that of the PSE and the reduction seems to be 
ective when the SNR is low. 

formance of the improved WVD 

brmance of the proposed method is illustrated for both FSK and chirp signds. 
red by the WVD, the signals are converted to analytic signals by computing the 
ransform of the original real signal. The Hilbert transform has been realised by 











Figure 3. Comparison of the resolution and ripple reduction by PWVD and MGDCM. (a) 
WVD slice by: - rectangular window (RW), — Hamming, - - MGDCM, (b) WVD slice by: 
- RW, — Hamming, -- MGDCM, (c) WVD slice by: - RW and -- MGDCM, using |iV(w)|2. 


For the FSK signal, the number of lags considered is 31. To reduce the crossterms, 
i-point boxcar smoothing is applied along the time axis, for the autocorrelation function £ 
;ach lag. Further, a discrete FT of length 128 points is used. For the PWVD, the SIACR i 
veighted by a Hamming window. For the proposed method, to avoid the negative spectrt 
values in the PSD, SIACR at zeroth lag has been lifted by a factor of 100. The | A((i))| 
jstimate is obtained by considering the inititil 8 cepstral coefficients in computing th 
moothed PSD. 

To get a comparative performance view, a section of the TFR obtained by the propose 
nethod and by the PWVD, corresponding to the stune inst£mt of time, £ire shown in figure: 
rtiis brings out the ripple reduction ability and the frequency resolution preserving abilit 
)f the proposed method (figure 3a). Further, in the MGDCM domain, a better improvemet 
n frequency resolution is observed (figure 3b) and this is due to the very nature of the GE 
rhe use of the |V((w)p esti m ate is unable to remove the ripple (figure 3c) and this is du 
0 the fact that the fluctuations in it are different from what is required to be canceled i 
he MGDC (for the ripple on the floor). ' 

The TFR obtained by the proposed method removes ripples/ringing due to abrupt trur 
:ation and at the same time preserves the resolution of the rectangular window (figure 4b 
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Figure 6. TFR representation for FSK signal with white noise by (a) PWV 
and (c) MGDCM after second smoothing, (d), (e) and (f) Contour plots ft 
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they support the fact that the resolution is improved over the PWVD. If 
peak location is of interest, then the GDMCM, above a certain thresl: 
ted in the TF plane rather than the PSD. This provides additional free 
(figure 3b). 

For the linear chirp signal, the results obtained are shown in figure 
smoothed along the time axis by a 9-point boxcar window. For this cas 
lag is increased by a factor 100. A floor level above a minimum valu 
the performance. However, increasing the floor level by a very large 
feet the TFR estimation. For the estimation of | A(n;)|“, the first 8 cej 
are used in computing the smoothed PSD. Figure 5c shows the TFR 
proposed method with a second time smoothing (box car of 5 pointJ 
effect was less when compared to that of PWVD even prior to the s 
(figure 5b). However, the second smoothing helped in reducing the ric 
which occurs at the region of crossing of the two chirps (figure 5c). Th 
crossterms may be due to their nature and as they are similar to ripp 
quency axis, they are flattened. In this case also, the frequency resoli 
and the ripple effect is reduced. The contour plots shown in figure 50 
results. 

The performance of the proposed method, for the above signals in the 
at SNR = 5 dB, is shown in figures 6 and 7 respectively. In both cases, 
method, the spectral peaks due to noise are reduced relatively to those oi 
is expected as the MGDCM not only removes the zeros due to the ripj 
the zeros due to noise and it does not distinguish between the two. 


5. Conclusions 

A method was proposed for estimating the power spectrum for a comp 
combines the concepts of the GD for a complex signal and the modi 
this approach was applied to the WVD to remove the ripple effect, ( 
truncation of the autocorrelation sequence, by tailoring the modified C 
signal, MGDCM. 

The proposed MGDCM method provides a spectral estimate which 1 
lower variance than that of the periodogram, without compromising on S] 
Depending upon the signal scenario, the proposed estimator has a varia 
about 70% to 2% of that of the periodogram. Further, as the variance 
effective when the SNR is low, its immunity to noise is very significant 
The improved WVD based on the MGDCM was found to be very effe 
the ripple effect due to truncation and in preserving the resolution of a re^ 
as no common window function is used. Further, since in addition to re 

due to the sneetml rinnle the yernc due tn u/hite nnice arp* jilcn rpmnvpd i 
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Abstract. A diverse choice of biquad designs is available for adaptive note 
filtering or line enhancement (ANF/ALE) applications. We consider one sue 
biquad proposed over a decade ago by David and coworkers and show il 
equivalence to two other formulations in terms of their pole-zero pair locatior 
in the z-plane. By equivalence, we imply that all these adaptive HR filtei 
possess comparable asymptotic performance. Further, a simple modificatioi 
involving only a normalizing factor, is seen to enhance the performance of th: 
filter. This modified ANF/ALE design is then shown to be equivalent to tw 
other recently proposed adaptive biquads. 

Keywords. Adaptive HR filters; adaptive line enhancement; adaptive sign; 
detection; adaptive signal processing; HR digital filters; poles and zeros. 


L. Introduction 

Adaptive HR notch filters (ANFs) provide superior performance at a lower computi 
ional cost relative to their FIR counterparts in suppressing tonal interference in wide 
>and signals. This advantage holds equally well for adaptive line enhancers (ALEs) i 
veil, where a sinusoid immersed in white noise needs to be extracted. As demonstrate 
n figure 1, the twin tasks of notch filtering and line enhancement may be viewed as con 
>lementary operations using a common framework. Further details of the topic can be ha 



Figure 1. Adaptive notch filter (ANF) block diagram. The structure can be ■ 
as an ANF or as an adaptive line enhancer (ALE), depending on the actual ou 
For sinusoids in white noise, a delay Z) = 1 is sufficient. 


different IIR adaptive filters have identical pole-zero pair locations in 1 
to yield comparable asymptotic (infinite data) performance. However, 
and tracking properties can differ depending on the specific parameter b 
adaptation as well as on the adaptation algorithm employed. 

As a reference design, we consider the biquad proposed by David and c 
cr a/1983; Hush et al 1986), hereafter referred to as the DEESHA filter (u; 
members in the team). This design is discussed briefly in the following ; 
then show its equivalence with other filters developed by Cho et al (19' 
& Martin (1991). In the fourth section, we consider a modified DEES^ 
improved performance and examine its relation to two other biquads pn 
& Martin (1989) and Regalia (1990). 


2. “DEESHA” ANF/ALE 


The DEESHA notch filter is given by 






/ 


[I - pz ^ -F r^z 


where 0 r < 1, and —2r < p < 2r. The factor p = (1 -|- r^) cos(£t 
acquire the notch frequency (on while “r” controls the notch bandwidt 
positioning of poles inside the unit circle. The bounds on p given above 
may be violated to handle very low (near d.c.) or very high frequenc; 
sampling rate) sinusoids without any penalty (Krishna & Hiremath 1995; 
of DEESHA filter suggested by Farhang-Boroujeny (1994) is 

H (r) = 1 

" 1 - (1 -|-ro)/0oz“' +roz~'^’ 

where po = p/d+f^) = cos(cu/y)andro = For a fixed ro, DEESHA. 
bandwidth notch filter. That is, its bandwidth is independent of the note] 
sinusoid in white noise. As a result, cuiv converges to the unknown fre< 




Asymptotic equivalence of some adaptive biquad notch filters 


Td2 ^3 



Figure 2. The dependence of pole locati 
on ro for DEESHA ANR The notch frequen 
is at = 7r/4 and the arrows indicate i 
creasing ro. 
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as a larger value of ro is selected within its upper limit, but at the same tiro 
filter convergence and tracking rates will also decrease. 

The notch filter requires its zeros to be located on the unit circle at an angle 
to the frequency of input sinusoid. The position of its poles is however detei 
the notch frequency as well as by the choice of kq. Figure 2 shows some pc 
ro is varied for notching at frequency 0.25 (normalized). In figure 3, the ani 
between a zero and its neighbouring pole is illustrated. Except for the cas 
with a normalized frequency 0.5, the poles lie at an angle closer to the ret 
zeros. At very low or very high frequencies (depending on ro), the poles ma}* 
real axis itself. This characteristic ensures that the biquad provides a consi 
notch irrespective of the input sinusoid frequency and thereby, a bias-free ( 
the presence of background white “noise”. 

The band pass filter (BPF) Hp, i (z) of the ALE corresponding to DEESl 
Hn, 1 (z) is given by 

rr I . (1 -ro)(po -Z~^) 

^ 1 - (1+ro)poz“^+roz“^' 

Note that the ANF and ALE are related by — 1 — ///?,! (z). T 

poles at the same location inside the unit circle as ANF, but has only a 5 
lying outside the unit circle. The zero can lie on the unit circle itself, if (jojsf 
SNR enhancement factor for DEESHA ALE is given by (1 + ro)/(l — r^) 
Figure 4 illustrates the magnitude response of the notch filter and line t 
on DEESHA biquad for some values of tq. 


DEESHA ANF / ALE Response 
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11 


5. DEESHA filter’s relatives 


iVe shall now examine the relation between DEESHA and Cho-Choi-Lee notch filters 


;Cho et al 1990). This latter design was derived from Rao-Kung filter (Rao & Kung 1984) 


rsing simple approximations and its performance was demonstrated (Cho & Lee 1993) 
;o be superior to the original Rao-Kung design in terms of reduced frequency bias. The 
R.ao-Kung ANF is given by 


Hn, 2iz) = 


\ a\z ' + a 2 Z 


(4) 


1 + aaiz~^ + a-a2Z~^ ’ 
sphere —2 < ui < +2 and 0 < 02 < I control the notch frequency wn = cos“'(—ai/ 
( 2 .^ 02 )), while a called pole contraction factor, places poles inside the unit circle on a 
radial line below the zeros for all notch frequencies. Note that the zeros will lie on unit circle 
anly for aa = L when this ANF becomes identical to yet another notch filter proposed bj 
Nehorai (1985). Unlike the DEESHA filter, both these ANFs are not constant bandwidth 
Rlters (i.e. bandwidth now depends on notch frequency) and thus suffer convergence bias 
induced by white noise component of the input (Krishna & Hiremath 1995; Regalia 1995), 
Bias can be a debilitating factor in applications such as phase jitter cancellation in carriei 
recovery phase lock loops (PLLs) used in high speed modems (Cupo & Gittin 1989). Che 
et al (1990) considered the lattice equivalent of (4) given as 

1+ko(H-ki)z“‘ 

Hn.2(z) = — -TT":- Z-ZT-. -ZT’ (5. 


1 + ?o(l A-q\)z-^ -\-q\Z 

where^o(l+^i) = a^o(l+^i),^i = ork\,kQ, = a\/{\+a 2 )md.k\ = a 2 -Then, assuming 
a pole contraction factor a 1.0, they used the approximations < 5 ^ 0 ( 1 +<?i) ^ kQ(l-\-ak\] 
and q\ ^ ak \, and substituted <70 = ko and q\ =: ak\ in (5) to obtain 

rr 1 +^0(1 +^i)Z"‘ 

Hn3{z) = , . , : • (6. 


1 +ko(l +aki)z~' +ak\z~^' 

For k\ = \, (which leads to zeros on the unit circle as does the Nehorai filter), this systen 
function finally simplifies to the Cho-Choi-Lee notch filter, 

1 -I- 2koz~^ +z~^ 


HNAiz) = 


1 -l-ko(l +a)2”^ +0CZ 


- 2 ’ 


(7 


where 0 < a < 1, and — 1 < kq < I. 

Cho & Lee (1993) termed the design in (7) as a “lattice ANF” and compared it 
performance with the Rao-Kung and Nehorai filters which were designated as “direc 
form”. This loose choice of terminology creates an impression that the Cho-Lee filter per 
forms better (in terms of bias) somehow because it is a “lattice ANF’. Comparing (7) witl 

T’AXU’TjOU’ a -Pi r M/-«• "D 11 7 •«^ /3k 1 ja-** i F 1 r> 1 fW TW1 O F/al ^ 
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Adaptive 

Algorithm 


eN(n) 


Other 

inputs 


Figure 5. ANF/ALE based on Mukund- 
Martin’s Resonator in a loop (RIL) structure. 


inside a stable feedback loop having a loop gain (or, damping factor) ’G’. Mukund 8c 
Martin explored its use mainly for line enhancement and for estimating the unknown 
frequency cos of a sinusoid but it can be used for notch filtering as well. The resonator, 
which is a BPF with infinite Q-idiCiox is given by 




I 4 

z zr z z^ 


lajtz ' - 2z ^ 
1 - 


( 8 ) 


where zr = exp(jo)R) oraR = cos(cor) determines the resonant frequency cor. The filter 
parameter ( 2 ft is adapted (Mukund &. Martin 1991) to obtain cor — cos upon convergence. 

To ensure loop stability, we require 0 < G < 1/2. Using straightforward algebra we 
can derive the system function between the input xin) and notched output e^in) as 


fiN.siz) — 


1 


lORZ ‘+ z ^ 


1-2(1-G)aftz-' +(1 -2G)z- 


(9) 


Comparing with (2), this is identical to DEESHA ANF for G = (1 — ro)/2 and ur = po! 
In fact, the performance of Mukund 8l Martin’s “resonator in a loop” design depends on 
the selectable value of loop gain, ’G’, to the same extent as DEESHA ANF depends on 
ro- Note, however that the notch filter bandwidth reduces as G approaches zero. 

Similarly, we can also see that the input x{n) and resonator output xs{n) are related 
by a system function, z~^ Hf^siz) where Hp^siz) is identical to the BPF of DEESHA 
ALE given in (3). This equivalence between DEESHA filter and the more recent Mukund- 
Martin’s (1993) design may appear surprising or even disappointing, if it were to be naively 
expected that the infinite-0 resonator in the loop would offer arbitrarily high frequency 
selectivity. Mukund & Martin (1991) further showed that their ALE provides an SNR 
improvement factor of (1 - G)/G, which exactly corresponds to (1 -f ro)/(l — ro) given 
by DEESHA. 


4. Enhanced “DEESHA” ANF/ALE 

For moderate values of ro (usually required during initial convergence or tracking time vary¬ 
ing frequencies), the DEESHA filter performance is not satisfactory. The SNR improvement 






uj aunit: uuuyiivt^ uiCfUULi nuicriJlucres 


/y 

tor for ALE (1 + ro)/( 1 — ro), remains small. Also, the magnitude response of the ANF 
ay from its notch frequency exceeds unity (figure 4). In fact, the response at d.c. and 
lormalized folding frequency, / = 1, equals 2/(1 + ro). Both these problems can be 
iply solved using the structure shown in figure 6. Here, the new notch filter output 
E (n) is only a scaled version of the original output (the notch filter pole-zero locations 
;he z-plane remain undisturbed). The scale factor ’ S’ is selected such that the magnitude 
ponse at / = 0 or 1, equals unity. That is, we require 5 = (1 -h ro)/2. The modified 
ignitude response is plotted (figure 7) for two different values of ro and may be compared 
th the earlier responses. 

The scaled output e^Ein) of the notch filter is then used to derive an enhanced ALE 
;nal xsEin). The enhanced ALE system function HpEiz) given by 


Hfe(z) = 1 - Hne(.z) — 


(1 - ^o) _ 1 - z ^ _ 

2 1 - (1 - ro)poz“’ + roz"- 


( 10 ) 


3 its poles at the same location as Hpiz) in (3), but has two real zeros, that remain fixed 
z = H-1 and z = —1 on the unit circle. In particular, these fixed zeros are of great help 
ten ro is small (i.e. the filter bandwidth is large). This is brought out in figure 7 where 
; magnitude response of the enhanced DEESHA filter is plotted. It can be easily shown 
It the SNR enhancement obtained using (10) equals 2/(1 — ro). Depending on the value 
ro selected, the improvement over DEESHA ALE can range from 1 to 2 dB or more. 
We shall now consider two more ANF/ALE designs and draw a relationship with the 
hanced DEESHA filter. Both these designs also satisfy the relation HpU) = 1 — 
tween their BPF and notch filter as in (10) above. That is, there is no explicit delay ahead 
the BPF here, unlike the general framework of figure 1 (the “missing” delay, which 
essential, is actually incorporated within the notch filter transfer function). Finally, the 
^E components of these two designs offer an SNR improvement factor identical to that 
enhanced DEESHA. 

The first of the two designs was proposed by Kwan & Martin (1989). To conserve space, 
; look only at the BPF (in ALE) here, as the corresponding notch filter can be easily 


Input 



Figure 6. Modification to enhance DEESHA ANF/ALE performance. 




Enhanced DEESHA ANF / ALE Response 



Figure 7, Magnitude response of enhanced DEESHA ANF/ALE. The notch is at = 
0.25, and ro = 0.81 and 0.64. These characteristics may be compared with those in figure 4. 


derived: 


Hp^dz) = 


£2 _ (1 - Z _ 

2 [1 - (2 - C2 - cj)z~^ + (1 - C2)z~^] 


(H) 


Here, c\ tracks the frequency of input sinusoid and has to be adapted, while C 2 controls the 
filter bandwidth. A simple eyeballing of (10) and (11) suffices to bring out the relations 
''O = 1 — C2 and ro = (2 — C2 — c |)/(2 — C2). 

The other filter we consider here was proposed by Regalia (1990). Looking once again 
at only its BPF, we have 


Hpjiz) 


1 - sin(g2) _ (1 - z~‘^) _ 

2 1 - sin(ei)[l + sin(02)]z-^ + sin(02)z-2 ’ 


( 12 ) 


with obvious relations to enhanced DEESHA given by po = sin(0i) and ro = sin(02)- 
Q\ controls the resonant frequency of ALE (same as notch frequency copj of ANF) and 
can be adapted using any appropriate algorithm to acquire and track the input sinusoid. 62 
determines the filter bandwidth similar to the factor, ro in case of enhanced DEESHA filter. 


5. Conclusions 

In this correspondence, we have established some relationships between different ANF/ 
ALE designs. Looking back, while some of these inter-relationships seem obvious (in 
particular, between DEESHA and the “lattice filter” of Cho, Choi and Lee), surprisingly this 
















■ J 


~r 


5pect has not been highlighted elsewhere. Even in certain instances where two different 
Dtch filters have been examined, invariably the comparison is between a so called “lattice 
Iter” and an inherently biased design such as that of Nehorai’s (Cho & Lee 1993; Regalia 
995). 

Our curiosity to explore these relations was aroused by the pole-zero plots of these vari- 
us designs. We were assisted in this exploration by Farhang-Boroujeny’s (1994) recasting 
f the DEESHA filter function. As the links between the DEESHA filter and others be- 
ame clear, this naturally led to the simple modification of figure 6 which improves its 
erformance both in notch filtering and in line enhancement. 

What is the implication of these relations that have been shown between various biquads 
Dr ANF/ALE applications? Mainly we have the assurance that equivalent designs will 
ffer similar asymptotic performance, in terms of bias-free convergence or for ALE, the 
NR improvement factor. Also, we hope that the confusion compounded by the wide 
ariety of design choices is reduced by the equivalences discussed here. However, the 
onvergence and tracking behaviour can significantly differ depending on the adaptation 
ichnique employed and on the specific parameter being updated. For example. Regalia’s 
esign given by Mukund &l Martin (1993) where 6 \ is adapted, is quite suited to handle 
ery low or very high frequencies. As for the enhanced DEESHA biquad, either po or ruv 
an be adapted. Its acquisition and tracking behaviour for very low or high frequencies will 
e superior for the latter choice (Farhang-Boroujeny 1994; Krishna & Hiremath 1995) and 
/ill be on par with that of Regalia’s, All these issues require careful consideration while 
electing a particular biquad design and adaptation approach for a given application. 
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Abstract. The lattice filter has several desirable characteristics that make it 
attractive in adaptive applications. The Lattice Stochastic Gradient (also called 
Gradient Adaptive Lattice algorithm) is popularly used to adapt the lattice filter. 
Howerver, a theoretical study of the bias in the reflection coefficient and the 
convergence of the LSG algorithm have not been studied extensively yet. This 
paper presents some theoretical results on these issues. It is hoped that the results 
will also present some insights into the factors affecting the convergence of the 
filter. 

Keywords. Adaptive filters; lattice filters; convergence analysis; adaptive 
algorithms. 


Introduction 

le lattice filter is useful in adaptive applications because it possesses several desirable 
aracteristics such as orthogonalization of the input, faster convergence, simple stability 
ecking criteria etc. The Lattice Stochastic Gradient (LSG), also called Gradient Adaptive 
mice (GAL) algorithm, js popularly used to adapt the (FIR) lattice filter as it is fast and 
IS a rea,sonably low computation. 

In adaptive filtering literature, there is a lack of rigorous theoretical analysis on the 
sues of bias in filter coefficient and the convergence of the algorithms. The analysis of 
e convergence of the LMS algorithm for the transversal filter is handled by making the 
idependence assumption’, which though handy, seems somewhat artificial. An analysis 
hich does not use this assumption is detailed by Douglas & Meng (1992) for a restricted 
nd of input. 
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Such rigorous analysis is absent in the case of the lattice filter. Honig & Messerschmitt 
have proposed convergence models by Honig & Messerschmitt (1981) for analysis of 
speed of convergence that are partly analytic. 

This paper attempts to prove absence of bias in reflection coefficient and the convergence 
of the LSG algorithm using more rigorous arguments. The notation used is that followed 
by Honig & Messerschmitt (1984). 

The order update equations of the lattice filter are: 

■ ef{T/n) = efiT/n - 1) - ■ eh(T ~ l/n - 1), 

et{T/n) = CbiT - l/n - 1) - • ef{T/n - 1), (1) 


where ej and e/, represent the forward and backward errors respectively, T is the time 
index, kl and k^ are forward and backward reflection coefficients of the nth stage. In the 
LSG algorithm, kl = 

The lattice stochastic gradient algorithms is: 


k„{T + \)=kniT) + 


1 


E(J/n - 1) 

\ef{Tjn - l).e/,(T/n) + eb{T - Xjn ~ 1) • ef{T/n)], 

E{,T/n - 1) = A • £(r - l/n - 1) + e}(J/n -1)4- ejiT - l/n ~ 1), (2) 
with ef(T /O = eh(T /O) = u(T) and E(0/n — 1) = 4-ve constant Vn = I,..., N. 


1. Theoretical results on the LSG algorithm 

Both the results presented below require assumption of ergodic input u{T) to the lattice 
filter Middleton (1960) for the definition of ergodic process used). The proof on conver¬ 
gence of the LSG algorithm requires, in addition, the assumption of ‘independent enough’ 
input (defined later). 


2.1 Absence of bias in kn {T) in the LSG algorithm 


Theorem 1. The lattice filter coefficients kn{T) are asymptotically unbiased when the 
LSG adaptive algorithm is used, assuming an ergodic input, ie. 


lim £[fc„+,(r4-l)] = ^"^V 


Proof. The LSG algorithm implements the following expression recursively 


^ • Cbji - l/n) 

Tl=Q{ef(i/n) +el(i - l/n)} 


(3) 


Denote the numerator and denominator of the above fraction as NUM and DEN respec¬ 
tively, for convenience. 

Here, we have considered the case A = 1 only, which is usually taken for stationary 
input. Now define 


No = E[2-ef{T/n)-eb(T ~l/n)], 





Lit LI 11 KC- LVLLLII^^ JLLLL^/J 


tJw/ 


DQ = E[e}(T/n) + el{T-\/n)l (4) 

We can write j NUM = A^'o + n{T) and ^ DEN = Dq + d(T). 

If u{T) is assumed to be ergodic, and we assume that the coefficients kj, j = , n 

ave converged (to constant values), then ef{T / n)-ei,{T —\ / n) and e^-{T / n)+ejj{T — \ / n) 
an be shown to be ergodic processes. Then limr^oo jNUM = Nq with probability 1, 
nd limr-^oo jDEN = Do with probability 1, 
which means 


lim n(T) = 0, with probability 1, and 

T-^oo 

lim diT) = 0, with probability 1, 

T-^oq 

/here ‘probablity 1’ means that the probability of the event that lim'f_^oo«(^j. T) ^ 0 
5 zero {n{Qj, T) denotes the jth sample function of the process n(T)). Now 


fc„+(r + i) = 


NUM/T 

DEN/T 


Do 


1 + 


n(T) 1 
No 


1 + 


d(T) 

Do J 


(5) 


E[kn+l(.T + l)] = k^/_^^-E 


1 + 


n(T) 

No 


1 + 


cl(T) 
Do J 


( 6 ) 


/here k°^^ is the optimum Wiener solution in the MMSE sense. 

Now, we prove that for the expectation on the RHS, denoted by £[(•)] for convenience, 
im 7 _,.oo £[(•)] = 1 under the ergodicity assumptions made above. 


£[(•)] = £ 


1 + 


n{T) 

No 


1 + 


d{T) 

Do J 


jeS 


, , rt(Q;.T)-l 
No 

1 , d(Qj,T) 

* + Oo J 


R(%), 


(7) 


!£[(•)]-11 = 



rn(Qj,T) d{Q.j,T)-] 

E 

jeSc 

Nq Do 

, , didj/n 

L ^ Do J 


• Pi^j) 


+ E 

jeSc 


rn(n;,r) 


No Do 

• Pi^j) 

, , di0.j,T) 

L ^ J 


sE 

jeSc 


n(Q^j,T) _ d{Slj,T) 

A'o _ Dp 

, , d(Qi,r) 

^ Do 


■Pi^j) 


+ E 


n(Qj,T) _ d(,Qi,T) 
Nq _ Dp 


1 + 


d{aj,T) 

Do 


■ Pi^j) 


vhere S denotes the joint sample space of the random processes niT) anddiT), Sc denotes 
he subset comprising all sample functions , T) and d(f2y, T) which converge to zero 
IS T 00 , and Sc denotes the complement of Sc- 
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Considering the second summation in above inequality, we see that since [ 1+d , T) / 
Dq] = DEN/Dor, and DEN is always positive non-zero, hence we can bound the absolute 
value of the fraction by an upper bound M Then, 

jeSc j(Sc 

which is zero under the ergodicity assumptions and above. 

Considering the first summation in same inequality, we see that since all n(Qj, T) and 
all d{Qj, T)eSc converge to zero as T -> oo, hence, given a 5 > 0, we can find a 7b such 
that |n(^', T)\ and |<7(%, T)\ ar both < 6 for 7 > Tq. 

Therefore, for 7 > 7o, the absolute value of the fraction in the summation < k-S, where 
k has a finite value. Hence, 

J2[(-)]-P(Qj) <k-S- ’ 

jeSc jeSc 

I by ergodicity 

which implies that 

\E[i-)]-^ <k-8, for7>7o. (8) 

Hence, limr^oo DfC-)] = 1 which implies that 

lim E[kn+^iT+\)]=kl;’:l^. (9) 

r-s-oo 

By following a method similar to the one used above, if we consider 


£[{fc„+i (7-1-1)-/:,71',}'] 



^( 7 ) \ 
Do ) 



/ = 2, 3, 


then, it can be shown that the term £*[(•)] on the RHS converges to zero asT oo. 

Thus, all central moments of kniT) 0 as T oo, which means that kn (T) converges 

to a constant value. 

Thus, under the assumption of ergodic u (T), it is proved that k] (T) converges to a fixed 
value which is also the optimum value i.e. k\ obtained has no bias. 

Now, ergodic w(T') and convergence of k\{T) to optimum value is sufficient to prove 
the convergence of k 2 {T) to its optimum value, and so on. 

Thus, all kn{T) converge to their respective optimum values. □ 


2.2 Proof of convergence ofk^iT) to a stable stationary point when the LSG algorithm 
is used 


Theorem!. Assuming an ergodic and independent enough* input (defined below), the 
lattice filter coefficients kn {T) converge to stable stationary points, when the LSG adaptive 
algorithm is used. 


roof. This proof uses the method of conversion of a recursive difference equation to an 
rdinary differential equation (see Goodwin & Payne 1977, pp 192-196). 

Consider the LSG algorithm (2). It can be rewritten as 

kn{T + \)=kn{T) + —-- . [2 • ef{T/n - 1) 

£(r/n — 1) 

■eb{T -\ln-\)-kn{T) 

.{e}{T/n-\) + el{T (10) 

T 

E{T/n - = - l) + 40' - I/« - !)}• (11) 

i=o 

Assume ergodic input M (7). Also, assume that the stages 1,..., n —1 have ‘adapted’; by 
'hich we mean that [ki (T)], i = 1,..., n — I have each converged to a stable stationary 
oint. 

We prove that in such a case k,iiT) will also converge to a stationary point 

For convenience, denote expectation operator £'[(•)] as (•)• 

From (10), we can write 


M-\ 

kn{T + M)^kn{T)+ Y. 

7=0 


1 

.£(r + i/n-l). 


•[2 • ef{T + j/n - 1) • eb(T + y - 1/n - 1) - /t„(r + j) 

.[ej{T + j/n - \) + eliT + j - l/n - m- ( 12 ) 


et 7 = 0 denote the time when all the first (« — 1) stages have ‘adapted’. Then from 
[ 1), we get 


'here 


_ ^ 

£(7/n-l) = X^l^-^'.{2.P„_i}, 

7=0 


Pn-\ = e^fUIn - 1) = e'^bU - l/« - 1)> and 7„_i 

1 independent of time by the assumption that fc/(y), i = 0,..., n — 1 have converged to 
able stationary solution. 


=^EiT/n-l) = 
^EiT/n-\) « 


1 - A 
ln(l/A) 


2-Pn 


1 - A 


■ 2 • Pn—l 


7. 


constant=:j6 


(13) 


he approximation holds when 7 <§( 11/ ln(A)|. 

Now,since£'(7/n —1) = /S-7, we use the heuristic that for large 7, E{T/n — \) = f-T 
Iso. Then, since increment in kn (7) is inversely proportional to (7 + y). for 7 » M, we 
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can replace + y) in (12) by kniT), which gives 


M-\ 

kniT + M) = k„{T)+ Y, 
7=0 


1 

.EiT + j/n-l). 


denote by ?RO{T+j) 

■ [2 • e/iT + j/n - 1 ) • eh(.T + j - l/n - 1 ) 

-kn (T) - {ej(T + j/n - 1 ) + eljT + j- l/n - 1 )}] 

denote by SUM(T+y) 


M~\ 

k„(T + M) = kniT)+ 

7=0 


1 

,E{T + j/n-l), 


•[PRO(r H- j) - kn{T) • SUM(T + j)l 


( 14 ) 


Now, let us assume that the input u{T) is, in addition to being ergodic, also an ‘in¬ 
dependent enough’ process, by which we mean that 3 a finite ‘Z’ such that u{T) and 
u{T — j) are independent Vy > 1. We shall call this the assumption of ‘independent¬ 
enoughness’. 

Next, we note the input samples that appear in the various terms of the summation 
in (14) 

E(T + j/n — 1) has samples u{T + j), u{T + J — , u(0). 

PRO (T -I- J) has samples m(T -I- J),..., u{T + j — n) and ki, i < n — 1. 

SUM (T -f j) has samples u{T -f j),..., u{T + j — n) and ki, i < n — 1. 

kn(T) has samples u(0), m(1), ..., u(T — 1). However, by the heuristic that EiT/n — 
1) oc T, it seems reasonable to believe that the dependence of (r) on u (t) will ‘decrease’ 
as T increases. At this stage, we cannot rigorously justify this, and hence, this decreasing 
dependence can be considered an assumption. We call this the assumption of “decreasing 
dependence”. 

Since we are considering larger, we can write ECr -t- j/n — l) ^ — 

1) -i-e^(T — l/n — 1)], so that EiT + j/n~l) now has samples uiT -f- j —n — l),u(T + 
j - n — I - , u(0). 

Comparing the samples u(v) in E(T + j/n — 1), PRO(r -I- j) and SUM(r -t- j), and 
noting that fy, Z < n — 1 have already ‘adapted’ before T = 0, we conclude that by the 
assumption of ‘independent-enoughness’, E{T -\- j/n — 1) is independent of PRO(r -t- j) 
andofSUM(r + j). 

Again, comparing the samples m(t) in kn (T) and PRO(r + j), SUM(r -I- j), we can say 
that for j > (I + n - kniT) is independent of both.PRO(T -|- ;), andof SUMCE -I- j). 
Choose M large enough (but <?C 3"), such that (•) ^ Therefore, for 

all terms in the summation, kniT) is independent of PRO and SUM. 

Finally, if £(7 4- j/n — l) ^ — 1)-|-c|(t — l/n — 1)] where p is large 

enough so that the assumption of “decreasing dependence” of kn (T) is valid for terms of 
EiT + j/n — 1), then kniT) will also be independent of EiT -f j/n — 1). Under these 
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nditions, we can write (14) as 


w-i 1 f 

kn{T + M)=kn{T)+ ^T—T- 

/=0 ^ ^ '- 


1 


EiT + j/n-DliT + j) 


r). 


■[PRO(T + J) - k„(r) - SUM(r + ;•)]. 


Since u(T) is assumed to be ergodic, hence e-^ {Jjn — 1) + e^ij — 1/n — 1) is also 
;odic. Therefor e, according to lemma 1 (see appendix A), we see that for X = 1, 
[E{-)I{T + J)] l/[E(-)/(T + j)] as r —> oo. We assume that for large T (but 

|l/ln(A.)|), we can replace l/[E(-)/(T + j)] by I/{£{■)/(T + j)\ even for X ^ 


Finally, noting that PRO(7 + j) = ryf,, SUM(r + j) = 2 • P„_i and E{-)IT = yS, 
lere r/i,, Pn-i nd ^ are independent of time for the conditions of ergodicity of «(r), and 
laptation’ of the lower n — 1 stages, we can write the above equation as 


kniT + M) = kn(T) + ~ . [ryy, - k„(T) • 2 • P„_, 
P 


M-] r 


E 

7=0 


iT + j 


small 

mverting this difference equation to an ODE 

^(^«a)) = 7 • [r/h - knit) ■ 2 • P„-|]. (15) 

le ODE gives as sthe stationary point 


knit) = rfb/[2- P„_i] 

^ E[2 ■ efjT/n - 1) ■ ehiT - l/n - 1)] 
E[e}iT/n-l) + el{T-l/n-])] ’ 


lich is also the optimum Wiener solution in the MMSE sense. 
To check for stability of the stationary point, we write 


{Akn ( 0 ) — 

at 


^ 2 P . 

lE 


Afc„(f). 


;nce stability is guaranteed iff 2 - P„_]/y3 > 0, i.e. > 0 , which is true. 

Hence, we conclude that under the assumptions made, when the lower n — 1 stages of 
i lattice filter have converged to stable stationary points, then the nth stage will also do 
. The first stage does not require convergence of any other stage. 

Thus, the first stage ‘adapts’ the earliest. Next, the second stage ‘adapts’ (since the first 
ige has already adapted), and so on for higher stages. 

This completes the proof. □ 


Results and discussion 

bile the two results presented above can both be regarded as proving absence of bias in 
flection coefficient as well as convergence of the LSG algorithm, yet they are not quite 


identical. The first result holds only for the case k = 1 (which of course should be used 
for truly stationary input). 

The second result holds for other values of k also. The price to be paid is less rigour 
in the proof as several heuristics have to be introduced. However, we feel that the .second 
result provides insight into the working of the filter in the way that the first result does not, 
Foi*, example, a look at the ODE obtained shows that the speed of convergence depends on 
(1 — A.) / ln(l/A.), which is independent of the MSE Pn-] and is very nearly independent 
of k for values close to 1. 

4. Conclusion 

The paper has attempted to approach the issues of convergence of the LSG algorithm and 
the absence of bias in the reflection coefficient in a somewhat more rigorous manner than 
is available in literature at present. The conditions imposed on the input are admittedly 
very strict, but they would be satisfied by Gaussian processes, which are used in simulation 
studies. However, it is clear that mere wide-sense stationarity is not enough for convergence. 
It is to be seen whether the conditions imposed on the input can be weakened. 

The authors would like to thank one of our reviewers for bringing to our attention a 
recent paper by Fan & Liu (1993), which proves the convergence of the LSG algorithm 
under more general conditions. 

Appendix A 

In this appendix, we state and prove lemma 1. 

Lemma L lfw{T) is a stationary, ergodic random process with w^T) = wq, and E(T) = 
Tthen 

1 

U)o' 

Proof. The proof uses arguments similar to those in th proof given in 2.1 
Since w(T) is stationary ergodic, hence we can write, 

Em = wo + niT), 

where n{T) 0 with probability 1 as T —> oo. 

Hence, we can write 


Next, we show that the expectation on the RHS—^ 1 as T oo. Write 
1 _ n(T)/wQ 

l+n(7’)/u)oJ Li+«(T)/uio. 





lim - 

Too IE(T) 



lu^uiiiiunjur uuujJLivt; lunim jnitr.'i 


yi 


len, by writing the expectation as two siimnaations, in a manner similar to that used in 
5.1, we get 


1 - 

1 ] 

. 1 + n{T)/wo\ 



keSc 


ni^k< T)/wq 




I +n(Qk, T)/wq 
n{Q.k, T)/wo 


1 +n{Qk, T)/wq 


PiQk) 


■ P(^k), 


lere Sc consists of all sample functions niQ^T) of n{T) which converge to zero as 
00 , and Sc is the complement of Sc. 

Now, the second summation on the RHS is zero since the denominator of the fraction 
always positive an = 0- 

As for the first summation on the RHS, given any S, we can find a Tq such that for 
> To 


E 


nj^k, T)/wo 
I+n{Qk, T)/wo 


■ Pi^k) < P-<5 E 

A*6S<. 

1 


lere p is a finite number. 

Thus, we have shown that as T —co, the expectation on the RHS of (16) —> 1. 

This therefore proves the lemma. □ 
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imbalance correction in time and frequency domains 
ith application to pulse doppler radar 

N SIVANNARAYANA and K VEERABHADRA RAO 

SHL, Reseai'ch Centre Imarat, Vignyanakancha, Hyderabad 69, India 
e-mail: siva@rci.emet.in 

Abstract. Digital ln“phase(I) and Quadrature phase(Q) imbalance correc¬ 
tion schemes are presented for improving the balance between I & Q signals by 
rejecting the image frequencies due to imbalance. The imbalance en'ors in the 
analog and digital demodulation schemes are highlighted. Simplified coiTection 
schemes are presented for time and frequency domain imbalances. These cor¬ 
rection schemes are useful in radar and communication systems. In this paper, 
a digital I-Q scheme is also presented for a pulse Doppler radar, along with 
hardware configuration for implementation. 

Keywords. Digital I-Q; quadrature demodulation; pulse doppler radar; time- 
frequency analysis. 


Introduction 

i coherent receiver, In-phase and Quadrature (TQ) phase signals are derived by demod- 
ting the Intermediate Frequency (IF) signal. The I-Q signals should match in gain and 
3hase by 90 degrees. In the analog schemes, IF signal is demodulated using in-phase 
I quadrature phase carrier (Goldman 1986; Liu 1989; Tsui 1995) and is sampled in two 
innels. In the digital schemes (Liu et al 1989, Tsui). IF signal is sampled and demodu- 
td using sampled cosine and sine of the carrier. These I-Q schemes ai*e explained in the 
rature towards simplification in the implementation (Levanon 1988; Liu et al 1989). 
wever, these schemes tend to develop amplitude and phase imbalances. In this work, 
are proposing simplified schemes for correcting the imbalances in time and frequency 
nains with application to pulse Doppler radar. 

’n radar and communication systems, signals are sampled for digital processing. As per 
quist criterion a signal must be sampled at a rate more than twice the bandwidth. In 
nplex (I-Q) sampling, the signals can be sampled minimum at the rate of bandwidth. 
I signals avoid blind phases in sampling. This complex sampling also improves signal 
noise ratio (SNR) by 3-dB compared to the only real (I-channel) signal processing 
wanon 1988). 
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The analog I-Q detector has a limitation of matching the gain & phase in the I and Q 
channels which can be achieved only upto a certain level (Goldman 1986). The mismatch 
creates an unwanted image frequency with an amplitude 24 dB down to the main signal. 
To achieve further rejection, the penalty in terms of cost is very high (Goldman 1986). 
Correction algorithms are incoiporated in the signal processors to reduce the imbalance. 

In the literature (Liu etal 1989; Tsui 1995), digital LQ schemes ai*e presented. Analog 
LQ imbalance corrections are explained by Churchill et al (1981) and Levanon (1988) 
and the calibration procedure is presented by Pien*e & Fuhrmann (1995). Hilbeit trans- 
forni techniques ai*e described by Oppenheim & Schafer (1975). Simplification of digital 
sampling is available by Brown (1979), Considine (1983), and Frerking (1994). Recent 
literature (Liu et al 1989; Tsui 1995) on I-Q, stresses the imbalance free demodulation. 
However, such demodulation creates an additional imbalance due to delay between I-Q 
samples for Doppler-shifted signals. For example, in radars, the image frequencies have 
to be rejected up to 50 dB to 60 dB down in severe clutter to signal levels of 30 to 40 dB. 
From the clutter signal, the image frequencies can be calculated, but ignoring of these 
image frequencies creates additional blind zones in the Doppler plane. So to avoid blind 
zones and false target detections, the image frequency level should be brought down to the 
noise background. 


2. Modeling I-Q imbalances 

I-Q channels are modeled in four ways depending on the presence of eiTor terms either 
in I or Q channels. The most commonly used signal with eiTors is (Churchill et al 1981; 
Levanon 1988) 

x{t) — a{t)[{\ + s) cos(cot) -f J smicot + 0)], (1) 

where a{t) is the envelope, 'co ’ is the angular frequency, 's' is the amplitude imbalance, 
'6' is the phase imbalance between I-Q channels and V’ is the time. The imbalance ratio 
is derived for (1), For simplicity a{t) is taken as amplitude A. 

x(?) = A[(l -)-£^) cos(<;uO "h j sin(cZ)r -t* 0)] 

= A/2 [(1 + sKej^' + 

= A/2 + e + + s - 

The image frequency strength X(—co) to main component strength X(co) is given by 

X{-co) _ (l+e-e-j^) 

Xio)) ~ i\+e + eJ^) 

[1 + £ — cos(0) + j sin(0)] 

[1 + e + cos(0) + j sin(0)] 

_ [s^ + 2s + j2{\ 4- £) sin((9)] 

[2(1 + £ + £2/2 + £ cos(0) + cos(0))] 

^e/2 + jO/2, for small values of £ and 0. 



imoaiance correciion in nine ana frequency aomains 


y:) 

ver ratio P-qjI is derived by taking |X(-rt»)p/|X((X))|^ and is given by (Churchill 
'/1981; Levanon 1988; Liu etal 1989) 

P-co/Pco^{e^ + 9^)/4. 

lilarly we can derive expressions for the remaining three cases. For the four possible 
)alance cases, the image to main signal strength is given by 


Signal model x(0 = 

X(-ft))/X(ft)) ^ 


A[( 1 -fa) cos(ft)t ) ~\~ j sin (ft)/ -1- (^)], 

(£/2 + ;0/2), 

(2a) 

A [cos (ft)?)-1- y(l+e)sin(ft)/-t-0)]. 

(-e/2 + y0/2), 

(2b) 

A[cos(ft)H-0)-F j(l-t-e) sin(ft)t)]. 

(-e/2-y0/2). 

(2c) 

A[(l “F a) cos(ft)/ -F 0) -F./ sin(ft)/)j. 

(a/2 - ;0/2). 

(2d) 


jalance eiTors can be estimated and corrected using frequency domain data while per- 
ning calibration. In order to calibrate, reference input is given at different frequen- 
! preferably covering the whole band. Depending on the selected model, the ratio 
-cu) / X (ft)) is considered with proper polarity. This is useful in automated calibration 
eceivers (Pierre & Fuhrmann 1995). 

Imbalance correction to the analog I-Q demodulators 

he analog I-Q, imbalances exist while demodulating the input using sine and cosine 
riers at intermediate frequency. Imbalance correction for analog I-Q detector is pre¬ 
ted below. 

Imbalance correction in the time domain 
i I-Q signals with imbalances of case-1, (2a) is modeled as 
/(t) = A (1 -I- e) cos(ftjr), 

Q{t) = A sin(ft)r + 9). (3) 

s can be corrected digitally by transforming the data as (Levanon 1988) 

I\{t) = EI{t), 

Qm = PI{t) + Q{t), ( 4 ) 

ere E = cos(6)/(l -t- s) and P — - sin(0)/(l -f e). 
kfter transformation, 

/l(r) = A cos(6) cos(ft)r), 

Q.\{t) = A cos(0) sin(ft>0. 

5 scheme mentioned above requires correction both in I & Q channel paths, 
ernatively, we propose to correct this imbalance in the Q-channel path only by 
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modifying (4), This modified scheme is suitable for easy implementation in the signal 
processors. 


= I (0 (remains same as input), 

Q\{t)^{P/E)I{t) + {\/E)Q{t), (5) 

3.2 Imbalance correction in the frequency domain 

In this scheme, analog baseband I-Q signals are sampled and the error correction is 
applied in the frequency domain. Normally, frequency transformation like discrete fourier 
transform (DFT) will be performed on the sampled data. The image frequency component 
appears after transformation. Hence, the correction is applied uniformly for all the bins. In 
the frequency domain, the frequency bins beyond fs /2 are complex-conjugates of the bins 
upto fs/2 and vice versa (Oppenheim & Schafer 1975) (/^ is the sampling frequency). 
Hence, both amplitude scaling and phase rotation, i.e, {e/2 H- j9/2) is provided uniformly 
on all the bins after conjugation operation. 

For 7/-point DFT, let X (k) denote the transformed output and X'^(k) its conjugate; index 
k goes from 0 to A/ — 1. The comection data XM is generated as 

XM(k) = (s/2 -f j9/2)X^(k)\ = 1, 2,..., /V - I, (6) 

and the corrected data XC(k) is given by 

XC(k) = X{k) - XM(N -k)\ * = 1, 2,..., - 1. (7) 

XC(0) = X(0) in all XC computation equations. 

With this type of correction, the image frequency component strength can be reduced 
by 30 to 40 dB below the main component strength. 


4. Digital I-Q generation 

In a digital I-Q scheme, the carrier (IF) is offset at half the bandwidth B/2 (Brown 1979; 
Considine 1983; Goldman 1986: Tsui 1995) and is sampled at 2B frequency. Demodulation 
is performed using cos(n7r/2) & sin(?t7r/2) which takes { — 1,0,1} for integer values of ‘n’. 
Thus, demodulation is a simple operation of changing the polarity of the input depending 
on { — 1,1}. In this scheme, cos{nn/2) is ‘1’ when sin(n7r/2) is ‘0’ and vice versa. So, Q- 
channel is a delayed version by a sampling interval. This delay creates an image frequency 
component. A fixed delay creates linear phase shift with frequency (Oppenheim & Schafer 
1975). In analog I-Q scheme, phase shift is constant. 

4.1 Imbalance correction in the frequency domain for digital 1-Q scheme 

The sampled signals when analyzed in the frequency domain, give linearly varying phase 
shift with bin number. Hence, every correction scheme requires the frequency bin depen¬ 
dent phase correction. For a fixed sampling clock, correction coefficients are generated as 
function of the bin number. 



Let fs be the sampling frequency for I-Q signals and delay between I-Q samples is 
Sampling interval 4 and St need not be related. Normally, St is less than tj. In pulse 
)ppler radars L is the pulse repetitive interval and St is equal to 1/2B for the signal 
npled at 2B. 

x(t) = A[cos(a)t) 4 - j sin(a;r + 

= A/2[iej‘^‘ + _ g-j(o)t+coSi)-^^ 

= A/2[ej^‘il + ( 8 ) 

le image strength X {—co) to main component strength X (co) is given by 

X{-co)/X(co) = (1 - e"-^'“^')/(l + 

= y [sin(a;30/(l +cos(a>3f))] (9) 

jQ)Stl2. 

St' is denoted as in the discrete frequency domain. Phase error for a bin of one 
3- fs/N) of -point DPT is given by 

<t>b=2Tt{fJN)St. (10) 

lase error for other bins is given by 

= k=0,\....,N 

XMik)=^j[smi<i>ik))/(l + cos(<I>{k)))]X*{k) (11) 

^jmk)/2]X*ik). (12) 

[uation (11) does not impose any constraint, and (12) limits the rejection for higher 
lays (St) between I-Q signals. Image-free components XC(k) are generated using (7). 
Windowing is applied on the time domain data to reduce the sidelobes. Figure la shows 
; input spectrum and figure lb is the XM(k). Figure Ic is the frequency reversal of XM 
d this is subtracted from input and indicated in figure Id, which indicates the image- 
rrected spectrum. The main component is at bin 60 and image is at bin 4. Blackman 
ndow is used for this. 

Due to windowing, frequency spread increases (Levanon 1988; Oppenheim & Schafer 
75). Let — 2 be the excitation bin and due to frequency domain broadening, N - 3 
N — I bins are also excited. For all these three bins, phase shift is (N — 2) <S>b. For 
1 - l)th bin, phase correction weight is (N — 2)^b- But, the correction applied is given 
(N — l)<l>i, in (12). The error rejection is ((N — \) — (N — 2))/{N — 2)^ (1/(V), 
nee, at higher bin numbers, the rejection improvement is limited to 20 log(77). 

For example, using a 64-point FFT, the possible rejection is around 20 log (64) = 36 dB. 
ior to correction, the image frequency ratio is of the order of —14 dB for the bin number 
I (image bin at 4) and — 18dB image for bin number 61 (image bin at 3) as shown 
figure la. After correction, the image is reduced by approximately 36 dB. With this 
:hnique, the image to main component strength is —54dB as shown in figure Id. All the 
aphs are normalized with respect to the input spectrum maxima. Smaller values of the 
ectrum are limited to — 80dB in displaying figure 1. 
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Figure 1. Frequency domain correction - magnitude plots. 
Phase is not indicated in the plots, (a) Input signal spectrum, 
X(k), main component at bin 60 & its image at bin 4. (b) Cor¬ 
rection spectrum, XM{k), this takes care of phase also, (c) Fre¬ 
quency reversed of XM{k). (d) Corrected spectrum. This is 
generated by subtracting XM(—k) from Xik), 


4.2 Imbalance correction in the time domain for the modified digital scheme 

The time delay between I-Q signals requires frequency dependent correction. As explained 
in the previous section in the frequency domain, coiTection vector is generated as fre¬ 
quency reversed conjugated spectrum with bin dependent coefficients. Thus, this forms 
multiplication in the frequency domain. Multiplication in the frequency domain reflects 
as convolution in the time domain. Convolution operation on the input gives correction 
signal in the time domain. 

The corrected spectrum is given by, 

- XC(k) = X(k) - XM(N -k); k = 1,2,,,,, N - 1. 

Let //(*) = y[sin(<i>(/:))/(l +cos(cI>(/:)))], 
and h(n)=^IDFT[H(k)l 

where IDFT denotes tlie inverse DFT operation. The number of significant coefficients in 
h(n) are much less compared to N, The image free signal jcc(«) is generated as, 

xc(n) = x(n) — jcm(n), 
xjn(n) = h(—n) O x*(n), 

where G denotes circular convolution. Circular convolution is a computationally involved 
operation compared to the frequency domain correction. Also, time domain correction 
scheme lacks the facility for auto-calibration of the receiver. 

While utilizing this scheme the following observations are to be noted 
(a) Filter consists of only few significant coefficients. Hence, the filter coefficients can be 
truncated to 5 to 8 coefficients, which saves computation time. 







imbalance correction in time and frequency domains 


99 



Figure 2. 


Analog I-Q scheme. 


Circular convolution is performed on the finite data block. 

The number of filter weights along with FFT size determine the possible imbalance 
reduction. 


I-Q scheme for pulse doppler radar 

mlse Doppler radars, canier at transmitted frequency is pulse modulated and received 
0 is first converted to the intermediate frequency(IF). While converting the signal from 
0 baseband, I-Q signals are generated as indicated in figure 2. In this scheme, carrier 
F is demodulated using and the output is filtered using a low-pass filter (LPF) 

tandwidth approximately equal to the reciprocal of the pulse width. I-Q signals are 
ipled using two Analog to Digital Converters (ADCs) at sampling interval less than or 
al to the transmitted pulse width (r). 

)igital equivalent of the scheme is indicated in figure 3. In this scheme IF is 
ited at /o(= 1/t) and sampled at 4/o. This signal is multiplied with a carrier 
[—y(2jr/on(l/4/o))] = cos{n7z/2) — js.m(n7i/2) and is suitably low-pass filtered 
down sampled to generate the required output sampling rate. 

^s an example, transmitter pulse width of400ns. is considered. The Steps involved here 
as below. 

Create IF (/o) at 2.5 MHz, one cycle is created in 400 ns. width. 

Sample at 10 MHz (4/o). 

Generate cos(n:7r/2) [10 — 1 0 1 0 — 1 0 ...] period of 4 samples. 

Generate sin(n;r/2) [010—1010—1 ...] period of 4 samples. 

Pass through LPF to pass signals up to -/o and remove components greater than 2/o. 
The filter should have maximum of 4 taps (equal to pulse width). 

Decimate the filtered output by 2 to 4, to generate the I-Q signals at the required 
sampling frequency. 



M 


Figure 3. Digital I-Q scheme. 
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5.1 Choice of filter weights 

Transmitted pulse is contained in approximately 4 samples. Matched filtering keeps the 
filter with samples comparable to the input. Filtering extends the signal due to convolution. 
Hence, the filter taps have to be limited to maximum of 4. Two simple filters selected are 
h = [I 1] or [0.5 I 0.5]. First filter is a Haar smoothing filter with minimum number of 
taps. Second filter is a convolution of [ 1 1 ] & [ 1 1 ] with normalization. The main advantage 
gained through this choice is the scope for multiplication-less implementation. 


5.2 Implementation 

Block configuration for digital I-Q is indicated in figure 3 and implementation in figure 4. 

In this, sampled signal at 4/o is demodulated using cos(n7r/2) & sin(n7r/2). This 
sampled carrier takes the values {-1-1, 0 & —1}. For multiplication by -f 1, 2’s comple¬ 
ment block passes on the input to its output. For —1, it is reversed in polarity using 
2’s complement operation. Multiplication by zero is performed in association with Regl 
by resetting Regl. Regl, Reg2, Reg3 outputs form 3-tap filter [0.5 1 0.5] (register is a 
latch/delay element denoted as ‘Reg’ in this paper). For multiplication by 0.5, adder 
inputs are given with a right shift. J2\ performs addition of first and the third sample. 

performs combined operation. This output is considered with suitable down sampling 
from 2 to 4. Decimation operation is indicated as a down arrow with M in the diagram. 
For the decimation factor M, PRI has to /tM/ 4, where I is an integer. With this condition 
for PRI an integer number of range bins are obtained and M = 4 coiresponds to one 
sample/pulse width. One sample per pulse widtli gives straddling loss and to avoid this 
loss, it is prefen'ed to keep M = 3. 

For a selected M, hardware can be further simplified by optimizing the number of 
registers and adders. The main advantage of this scheme is the feasibility of configuring 
the whole hardware after ADC in Erasable Programmable Logic Devices (EPLD) or Field 
Programmable Gate Array (FPGA) without any multiplication operation. Hence, this type 
of scheme reduces hardware complexity. 



Figure 4. Digital I~Q implementation. 
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] Imbalance correction 

1 explained in the previous sections, digital 1-Q requires frequency bin dependent cor- 
:tion. This has to be incoiporated in the signal processors. However, this imbalance 
rrection may not be required in the ground based and clutter Doppler locked pulse 
)ppler radars. 


1- Options 

) Sampled signal is demodulated using cos(/z7r/2) & sin(/2jr/2), passed through a low 
pass filter and decimated suitably. 

) Sampled signal is demodulated as above and decimated by 2 or 4 in synchronism with 
the demodulating waveform (Liu et al 1989). 

) Sample at highest possible frequency n(4fo) to reduce phase imbalance, where ‘n’ is 
greater than or equal to 1. 

) For pulse compression waveforms like linear FM, the block LPF can be merged with 
the pulse compression filter in digital pulse compression schemes. This I-Q scheme is 
hardware intensive and so the possible options ai*e to be considered depending on the 
application. 


Conclusions 

this paper, different types of TQ imbalance correction techniques are presented in 
ne and frequency domains. These imbalance correction schemes are easy to implement 
both hardware and software and avoid the use of matched LQ detectors which ai'e 
stly. It may be easier to implement these schemes in the frequency domain and the 
jquency domain data can be utilized for amplitude and phase imbalance estimation in the 
tomated calibration procedure. This is an added advantage. The computations required in 
s above schemes are negligible compared to computation involved in DFT/FFT and other 
jnal processing operations. For digital I-Q, it is prefeixed to go for frequency domain 
irrection. Digital I-Q schemes are useful for radar and communication applications, 
lese schemes with higher image rejection will be useful for air-borne Doppler radars, 
le optimum hardware can be configured for the desired parameters depending on specific 
plication. 


le authors would like to thank Prof. V U Reddy, for his invaluable suggestions. The 
ithors would also like to thank Dr K V S S Prasada Rao and their colleagues for their 
pport and encouragement. 
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terpolation of erasure bursts via cosine-modulated 
terbanks 
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Hyderabad 500 082, India 
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Abstract. A novel and low-complexity approach for reconsti'ucting periodic 
erasure bursts in data sampled at greater than the Nyquist rate, using cosine 
modulated filterbanks, is described. In the case of interpolation of erasure sin¬ 
glets or doublets periodically repeated over 2M samples, the cosine modulated 
filterbank approach is shown to have a lower complexity (for a given restora¬ 
tion error) than a standard FIR interpolator. In the case of erasure triplets or 
quadruplets, periodically repeated over IM samples, the restoration error is 
primarily related to whether the M-channel filterbank’s stopband suppression 
is better than the condition number of a 2 x 2 matrix, where M is determined 
by the oversampling factor of the data. While the method used for erasure 
triplets and quadruplets extends to arbittary erasure bursts, the condition num¬ 
bers of the associated (larger dimension) matrices deteriorate rapidly with 
the increase in erasure length, posing practical problems such as the design 
of very high-attenuation filterbanks and lai-ge required implementation word- 
lengths. 

Keywords. Chebyshev filters; discrete cosine transforms; equiripple filters; 
M-channel filterbanks; quadrature mirror filters; signal restoration; signal sam- 
pling/reconstmction. 


Introduction 

; restoration of erasure bursts in oversampled multimedia data has assumed increased 
nificance in the context of cell loss in ATM networks. Many heuristic methods, such 
linear interpolation, repetition, or muting, are recommended in many multimedia stan¬ 
ds. Other systematic approaches, however, have been considered in the recent years, 
die Vaidyanathan & Liu (1988) examine the feasibility of correcting long error bursts 
he context of the sampling theorems, Marks and Radbel (Marks 1983; Marks & Radbel 
14), by considering the condition number of sub-matrices of a DFT matrix, exposed the 

103 


lU'-h 


o juyusirnnu ana kj iiircrnuin 


limitation^ of the DFT method of signal restoration. The cosine modulated filterbank ap¬ 
proach to burst erasure restoration described herein shows similar limitations. However, a 
poor condition number of a submatrix of a DCT matrix may be overcome by the design of a 
high stopband attenuation filterbank (in terms of achieving the desired restoration error). It 
is also shown that the approach proposed applies even to very slightly oversampled signals 
by increasing the number of channels M, albeit at the expense of the condition number 
of a submatrix of a DCT matrix. The method described is limited by practical considera¬ 
tions such as the design of very high stopband attenuation filterbanks and large required 
arithmetic word-lengths (for longer erasure bursts in signals approaching the Nyquist band¬ 
width). The method may be used to correct erasures either for isolated bursts (separated by 
at least the order of the prototype filter) or for either M- or 2M-periodic erasure bursts. The 
method described is also computationally efficient and in the simplest single 2M-periodic 
erasure case, obtains a factor of 2 in reduced complexity as compared to an ordinary FIR 
interpolator. 


2. New restoration method 

The restoration method uses a typical M-channel maximally decimated filter bank as 
shown in figure 1, where Hk{z) and F^(z), 0 < k < M — 1 are analysis and synthesis 
filters, respectively. 

Here the cosine modulated filter bank is used which is attractive with respect to 
implementation cost and design ease. The cosine modulation may be either of type II 
or type IV^. The impulse response of the analysis and synthesis filters hif^in) and 
fk(n) are cosine modulated versions of the prototype filter h(n). For type IV cosine 
modulation: 

hkin) = 2/z(«)cos (yik + 1)^ - (Ik + 1)^^ , 

fkin) = 2 / 2 (n)cos ^(2* + 1)^ (n - + (Ik + 1)^^ , 

0 </i < iV- 1, 

' 0<k < M- 1, 


where N is the length of h(n ). Here the cosine modulation uses a phase shift of (2/: 4- l)7r /4 
and satisfies the alias cancellation constraint given (Vaidyanathan 1993). The analysis 


^Quoting from Strang & Nguyen (1996): “The Discrete Cosine Transform (DCT) improves on the DFT for the 
Same reason that symmetric extension improves upon periodic extension. The symmetric extension is continuousT 
^The kernel for the modulation matrix is 



Yn/ , ^ 


cos 


1 &k 


where % is the phase shift determined by the alias canceling constraint (Vaidyanathan 1993). For type II modulation, 
V’ is the length of the filter and for the type IV modulation, N' is the order of the filter. 
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Figure 1, M channel maximally decimated filter bank. 
r with these impulse responses can be separated into 2M polyphase components 


2M-1 W2M)~1 

E E ^h(2qM^i) 

i=0 q=0 

X cos + 1) (^i - --Y~) ^ j) 

2A/-1 

= GiC-z^") 

i-o (2) 

X cos ((2i + 1) (i - ^ - (2t + 1>|) s-', 

iN/2M)-l 

where ^ 2 ;z( 2 ^Af + z)(-z)- 29 W 

q=0 

0<q < {N/2M) - 1, 

0<k<M 

an efficient implementation of this filter bank, the length of the prototype filter N is 
imed to be an even multiple of M, i.e., N = 2mM, where m is an even integer; this 
dition is not restrictive as a prototype of any length can be padded with an appropriate 
iber of zeroes. This simplifies the filter transfer function to: 


0<k<M-l. (3) 



ast algorithm for implementing this filter bank using the IDCT-IV like modulation 
rix (the modulation matrix used here are same as the DCT-IV matrix without the 
al scaling factor for the first coefficient) can be obtained by reordering the polyphase 
iponents as follows: 

^ f z-^/2g,+(m/ 2), for / = 0.1,... (3M/2) - 1, 

' [-z3^/^G,_(3m/ 2), for / = 3M/2, ..., 2M - 1. 
















In terms of this new sequence, the analysis filters are expressed as follows: 
M-\ 

Hkiz)= 

i=0 

xcosf(2fc + l)(2/ + l)^V 0<^<M-l. 

V 4-M / 

In matrix notation, the analysis filter bank vector h(z) becomes; 
h(z) = Tg(z), 


where 


(5) 

( 6 ) 


Hz) = 


■ Hoiz) 
Hdz) 




_ _ ^ ^ 

g(^) = 

and T is an M X M IDCTIV like modulation matrix with elements: 

n 


t]^i ~ 2 COS I {'Ik H- l)(2z + 1) 


4M 


The coiTesponding synthesis filters can also be expressed in terms of the DCT IV 

like modulation matrix: 


M-\ 


Fkiz) = Y. 2z-^iz-‘G'i(-z-^^) - z 






i-zr^^)) 


/=0 


X cos ( {2k + 1)(2/ + 1) 


TT \ 


am) 


), 0 < /: < M — 1, where K = N — 2M, 

(7) 


In matrix notation, the synthesis filter vector f(z) is written as 

f^(z)=z-"g^(z-^)T^, (8) 

where, f^(z) = [FqC^) Fj (z) • • • Fm-i (z)]. 

Similarly, the analysis filter and synthesis filter vectors for type II modulation kernel 
can be expressed in terms of the DCT II modulation matrix as given below: 

h//(z) = T//g//(z), (9) 

where 

■ HoU) - 

miz) 

5 


h//(z) 
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iki = 2cos (^i2k+ • 

synthesis filter vector t(z) is written as 

f[/(z) = 2“^g//(z"')T[/ (10) 

Dlyphase structure, similar to the one given, that implements the analysis filterbank for 
: IV modulation is shown in figure 2. 



Figure 2. Cosine modulated analysis filter bank in terms of IDCTIV. 






The analysis filters H\^(z) channelize the input signal into M subband signals Sk, which 
are in turn decimated by M. These subband components can also be expressed in matrix 
form in terms of the IDCT matrix and polyphase components as: 


S = TP 


(H) 


where S = [^o, 5],..., is a vector of subband components and P = [Po< Pi, 

..., Pm-\]^ is the vector of polyphase components Pi (i.e., P = g(z)X(z), where X(z) = 
[x(z) z~*x(z) • • • z~^“'x(z)]^) and T is a M x M orthogonal IDCT IV matrix. For an 
oversampled signal with a nonnalized bandwidth of r {l.F/ ,, where F is the bandwidth 
of the signal and Fs is the sampling frequency), and the last K = \i\—r)-M — \'\ subband 
samples are not present. 

The p missing (or lost) samples are so aligned that a minimum number of polyphase 
components are unknown. In the case of a burst of erasures, the erased sample indices range 
from M — p/2 (mod M) to M + p/l (mod M), so that only \p/2'\ polyphase components 
are unknown. Then the number of channels of the filter bank M are chosen according to: 


M = 


' (P + 1) 
2 


(1 -r) 


( 12 ) 


where, [12 denotes rounding to an even number towards infinity. 

For the filterbank with DCT-II modulation mabix, the number of unknown polyphase 
components is [(p +1 )/21. Therefore, when the length of erasure burst is even the number 
of unknown polyphase components for the filter bank with DCT-IV modulation matrix is 
one less than that for the filter bank with DCT-II modulation matrix. 

Next, split the polyphase vector into a vector of known components Pj and a vector of 
unknown components P 2 , so that 


Si 



Tfil 

.p. 

-S2^ 


H 

_1 

- , 

H 

P 2 . 


(13) 


where, Tj’,, Tf 2 , T^j and T 22 are the submatrices of a permuted IDCT IV matiix. The 
columns of the modulation matrix T are permuted according to the known and unknown 
polyphase components. As the input signal is bandlimited, the subband component vector 
S 2 can be set to zero. This results in a system of linear equations which are solved to 
determine the unknown polyphase components P 2 . 


[Tf2]P2 = -[T^,]Pi. 


(14) 


As shown in figure 2, when the erased sample burst indices range from M—p/2 (mod M) to 
M + p/2 (mod M), the first p/1 components of the polyphase vector P 2 /r = [Po» Pi. • • •, 
Pp_i] will have to be determined for the first M-sample input block. When the next 
M-sample block is the input, the last p/2 components of the polyphase vector P 2 L = 
• • •. Pm~-\] are unknown. The unknown polyphase component vectors 
alternate in this manner until all the erased samples are restored. 

In the preceding paragraphs, we have considered 2M periodicity in determining M from 
the given normalized bandwidth of r and burst length p. An alternate view of figure 2 is 
to consider a periodicity of M in determining M from the oversampling factor and burst 
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;th. In this case, each missing polyphase component in the top M polyphase components 
has a counteipart, separated by M from it, that is also missing. In this case, 
DCT submatrix is always constant for every M-sampIe block during the restoration 
:ess. The number of subbands M is determined by 


M = 


z^ + r 
1 +r 2 ' 


(15) 


number of channels of the analysis filter bank obtained using (15) is less than twice the 
iber of channels obtained using (12). Therefore, the interval between the enor bursts 
g 2M-periodicity (with M given by (12)) can, in some cases, be reduced by using 
)eriodicity (with M given by (15)). 


Error in restoration of erased samples 

estimated polyphase components are passed through the polyphase structure of the 
hesis filter bank to determine the erased samples. Therefore, the error in the restored 
pies is given by: 

\\SX\\ = l|g.(2"‘^P2(z)ll < Ilg.(2-')1!1|5P2(Z)II (16) 


TC 


¥P2(z)\\=^ 


8V2F z ^SP2F + * 

• z 3 P 2 /r 

8P2L + z~'^8 V2 l + •' 

> • Z SP2L 


\\ge(z II is the norm of the polyphase components in g(z) corresponding to P 2 F and 


he error of the estimates of the unknown polyphase components using (14) is due to 
error introduced into (14) by the assumption that the subband vector S 2 is zero. The 
r in the restoration of P 2 , 8^2 satisfies: 

PP2ll</^(T^2)^lj!|^l|Plll U7) 

;re, K (T 22 ) is the condition number of matrix T 22 ^ and 1| ||‘* denotes the norm of matrix 
vector (Kreyzig 1993). For error vector 5T;^2 be small, both the condition number 
and the norm subband vector S 2 should be small. If the condition number of the 
rix T 22 can be controlled (through the choice of M), the norm of the subband vector S 2 
be reduced by using a filter bank with higher stop-band attenuation.^ It should also be 
:d that as the condition number of the matrix T 22 increases, the word-length required 


condition number of matrix A/c (A), ivS the ratio of largest singular value of A is the largest singular value of A, 
m of matrix A is the largest singular value of A, 

uming a perfectly bandlimited signal in [0, (M - K)Tr), the norm of the decimated signal in any high 
and (that will be set to zero) has an upper bound of M - K times the supremum of the norms of the signals 
i first M — K subbands multiplied by the supremum of stopband gain of the K high subband filters. This again 
3e small if the maximum stopband gain of the prototype is small. 


M 

DCT IV 

k(T^2) 

P2i. 

DCT IV 

P2/rDCT 

TI 

P 2 L 

DCT II 

8 

91.2789 

11.2485 

24.4259 

41.5830 

12 

225.6398 

27.2672 

56.7647 

108.2118 

24 

954.6982 

114.6825 

231.7052 

472.3874 

32 

1711.200 

205.4472 

413.3550 

850.5666 


Table 1(b). Condition numbers for DCT II and DCT IV siibmatrices T 22 for different K with M = 8. 
T 22 for different M with K =2. 

K 

/r(T''2) 


^ *■22 



Pri/r 

PlL 

P2fDCT 

P2L 


DCT IV 

DCT IV 

II 

DCT II 

2 

91.2789 

11.2485 

24.4259 

41.5830 

3 

667.0022 

40.2057 

133.4108 

190.1040 

4 

1059.4000 

58.7428 

260.1650 

236.5874 


for solving (14) also increases and imposes practical constraints on use of the usefulness of 
the algorithm. In the case of a large erasure burst, the reconstruction (psuedo-QMF, perfect, 
or near-peifect) properties of the prototype filter has only a negligible effect on restoration 
error.® Consequently, in the following examples used to illustrate lost-sample restoration, 
we have used the near-equiripple pseudo-QMF bank designed using the method described 
in by Jayasimha & Hiremath (1998). 

The condition number of the matrix T 22 deteriorates with increasing M and with the 
number of missing polyphase components (i.e., the length of the lost-sample burst). The 
condition number of the matrix T 22 used while estimating V 2 F and P 2 L vectors is tabulated 
in tables 1 a and b, respectively. As can be observed from the table, the maximum condition 
number is for the DCT-IV modulation matrix. Therefore, when the condition number is 
large, the filter bank with the DCT-II modulation matrix is used (in particular, for odd p). 

To compensate for the ill-conditioning of matrix 11^211 should be small and the 
prototype filter must be designed to have high stopband attenuation. A method to design 
high-order, large even M, filterbanks (where M is an even composite number) is discussed 
by Hiremath & Jayasimha (1997). 


4. Interpolation of periodic erasure singlets and doublets 

A block diagram for the efficient restoration of periodic, erasure singlets, compared with 
the restoration of periodic erasure singlets using a type II maximally-decimated M -channel 
(where M is divisible by 4) cosine modulated filterbank, is illustrated by an example in 
figure 3 for periodic erasure singlets with periodicity of / = 4 and M = 12. 


^In pseudo-QMF banks, the stopband attenuation; however, the restoration error is of the order of the stopband 
attenuation 
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e 2. Restoration ESR for erasure bursts of various lengths (M = 8). 


ire 

Actual 

samples 

Restored 

samples 

ESR upper 
bound 

ESR 

practical 

let 

3352 

3353 

0.0026 

0.00029 

)let 

3666, 3352 

3668, 3345 

0.004629 

0.001465 

et 

3666, 3352, 2855 

3668, 3345,2873 

0.006688 

0.001317 

Iruplet 

3806,3666,3352,2855 

3779,3651,3352,2871 

0.005827 

0.005057 

tuplet 

3806, 3666, 3352, 2855, 2248 

3779, 3651,3352, 2871,2229 

0.007171 

0.005477 


he complexity per interpolated point for the restoration of singlets with M — 12, / = 4 
IV = 128 is given by C = I/M (DCTfy + DCTj^y + 3 • adds + DCTfj) + W/M 
iiply adds 25 multiplications +35 additions. 

he usual FIR interpolator derived from a 4th band filter (Adams 1991) for the restoration 
h similar performance) of these periodic erasure singlets has 89 taps and uses 67 
tions and 34 multiplications. 

I a similar manner, doublets can be restored using the modulation matrix of DCT 
or which no FIR interpolator is proposed by Adams (1991). Note that the condition 
her «-(P.^ 2 ) is one and the error in the estimation depends upon the subband component 

l- 

Erasure burst restoration example 

le randomly selected position in audio data sampled at 44.1 kHz (an example from 
e Straits’), bandlimited to 13.23 kHz, is used to test the described restoration method. 
( z= S, N = 128 low-pass prototype with specifications: Passband edge normalized 
uency: 0.01042 Stopband edge normalized frequency: 0.11458 Stopband attenuation: 
dB was designed (Jayasimha & Hiremath 1998). This odd-length prototype padded 
a leading zero, and used with a DCT-II modulation matrix for fast implementation, has 
:onstruction error of —54dB. The average restoration error to signal ratio, a measure 



Figure 3. Restoration of singlet s(M = 12, / = 4). 
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of the restoration’s accuracy, is: 


ESR = 


Y^ixin) -i(/i))^ 

n 



1/2 


(18) 


The restoration ESR’s for isolated bursts of various lengths 
table 2. 


for M = 8 are given in 


6. Conclusion 

A novel, low-complexity cosine-modulated filterbank approach to short-burst erasure inter¬ 
polation has been described. Potential applications areas of short-erasure burst restoration 
are in removal of lightning-induced impulsive noise in VLF and LF communications (de¬ 
tected by short-duration saturation of amplifiers and ADCs) and the extension of the ADC 
dynamic range when the bandlimited signal being acquired is slightly oversampled. One 
example of the latter application is that saturated signals in short regions of some 16-bit 
compact disk recordings may be restored and played back on 18-bit DACs. 
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Abstract. The scenario is that a bulk data transfer is being performed over a 
TCP connection, from a host on a local area network (LAN) to a mobile host 
attached to the LAN by a radio link. In an earlier work we had assumed that 
packet losses in a TCP connection over a radio link are statistically independent. 
In this paper, we extend this analysis to a Rayleigh fading link, which we 
model by a two-state Markov model. The bulk throughputs of TCP-OldTahoe 
and TCP-Tahoe are compared with and without fading, for various average 
signal-to-noise ratios. We also study the performance with a link protocol on 
the wireless link, and study the effect of varying the link packet size, the number 
of link packet attempts, and the vehicle speed. For the parameters of the BSD 
UNIX implementation, over a 1.5 Mbps wireless link, we find that, with fading, 
a signal-to-noise ratio of at least 30 dB is required to get reasonable throughput 
with TCP Tahoe or OldTahoe; this corresponds to at least 100 times more power 
than is needed without fading. 

For fixed signal-to-noise ratio, as the vehicle speed varies there are roughly 3 
regions of performance; at very low speeds (pedestrian speeds) the throughput 
is very good; at low vehicular speeds the throughput deteriorates, and again 
becomes very good at higher vehicle speeds. The speeds corresponding to the 
various regions depend on the parameters of the link protocol. 

Keywords. Mobile internet; mobile computing; TCP modelling; TCP over 
fading channels. 


iitroduction 

letwork scenario (figure 1) is motivated by the many recent experimental studies of 
performance over wireless mobile links (Bakre & Badrinath 1995; Caceres & Iftode 
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Figure L A LAN host with a TCP connection to a mobile host. 


1995; Balakrishan et al 1997). We are interested in the data throughput, from the LAN 
host to the mobile terminal, that can be achieved by various versions of TCP, when the 
wireless link introduces random packet losses. In our earlier work on this problem (Kumar 
1996) (see also Mishra etal 1993, Lakshman & Madhow 1997) a simple loss model was 
assumed: TCP packets were transmitted in their entirety over the wireless channel, and each 
packet was lost independently of anything else with probability p (i.e., the Bernoulli loss 
model). The Bernoulli loss model would correspond to a nonfading channel with additive 
white Gaussian receiver noise. It is well known (Jakes 1974; Parsons 1992) that when the 
mobile and the transmitter are not in line-of-sight, the mobile receiver antenna observes a 
superposition of multiple reflected and diffracted signals that are out of phase with each 
other. The phase relationship between the interfering multipath signals is continuously 
changing as the vehicle moves. Destructive interference of these signals leads to “fading”. 
The fade durations depend on the velocity of the vehicle-and the radio carrier frequency. For 
high speed radio transmission (e.g., 1.5 Mbps), typical carrier frequencies (e.g., 900 MHz), 
and the usual vehicle speeds, the fade durations are comparable to the transmission times 
of TCP packets, or are multiples of transmission times of typical link packets. Thus the 
packet losses cannot be modelled as being independent of one other. 

In this paper, we model the losses using a two state Markov model. We assume that a 
packet succeeds with probability 1 in the Good state, and with probability 0 in the Bad 
state. The TCP protocol alternates between packet transmission periods and loss recovery 
periods. It was shown by (Kumar 1996) that the TCP packet transmission process can be 
modelled as a Markov RenewaLReward process (Wolff 1990), the Maikov renewal instants 
being the epochs at which the first packet loss occurs in a transmission period, and the 
“reward” corresponding to successful packet transmissions. With the simple fading model, 
this Markov renewal-reward analysis is easily adapted. We obtain numerical results for the 
same parameters as by Kumar (1996), 

In a previous work, Gilbert has used a two state Markov model to model bursty error 
rates in digital transmission links (Gilbert 1960). Wang and Moayeri (1995) have studied 
multistate Markov chain models for radio channels, and have, in particular, developed a 
finite state Markov model for a Rayleigh fading channel. Recently, Chaskar et al (1996) 
have analysed a wide area TCP connection in which the receiver end system is connected 
to the wireline network by a wireless mobile link. A Rayleigh fading model is assumed 
for the wireless link, and a two state Markov model is used. It is assumed that the link 
protocol on the wireless channel repeatedly retransmits the link packets until they suc¬ 
ceed. Since this implies an increase in the effective transmission time of TCP packets, 
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lain concern in the paper is the additional buffer requirement (at the wireline-to- 
sss router) as a result of wireless channel losses. It is argued that for a Markovian 
g model (implying exponentially distributed fade durations) TCP will yield reasonable 
ghputs if the buffer space grows only as fast as the logaiithm of the bandwidth-delay 
ict. 

our analysis, we assume that there are adequate buffers at the wireless router so that 
r overflow is not a concern. On the other hand, we permit the possibility of packet loss 
i wirelesj; link, either because there is no link protocol, or because the link protocol 
1 the number of retries. The effect of losses on the TCP throughput is studied. We 
de results for throughput performance of TCP-OldTahoe and Tahoe as the signal¬ 
ise ratio varies, and compare the situations with and without Rayleigh fading. For 
signal-to-noise ratio, we then study the variation of TCP throughput with the vehicle 
1 . Varying the vehicle speed effectively varies the duration of the fades and the good 
is on the radio link, and results in interesting behaviour of the TCP throughput as a 
ion of the vehicle speed. 

is paper is organised as follows. In § 2 we develop the two state Markov loss model. 
I we describe TCP’s window adaptation protocol. In § 4 we adapt the throughput 
sis developed earlier (Kumar 1996) to the Markov loss model. In § 5 some numerical 
s and their discussion are presented. 


t Markov model for Rayleigh fading 


Review of the Rayleigh fading model 


adio carrier is digitally modulated; for Pulse Amplitude Modulation (PAM), the 
Tiitted signal is written as 


x{t) = \/2 Re 


Ske-’^'‘ pit - 


lc=—oo 


i p(-) is the baseband pulse, fc is the carrier frequency, and (^/t, %) is the complex 
ilating sequence (see, for example, (Lee & Messerschmitt 1988). Analysis of the 
position of multipath signals, in the presence of receiver mobility, yields the following 
1 for the received signal (Rappaport 1996): 


y(t) = VS Re 


Y • pit - 


(I) 


[/:=-oo J 

2 I'it) is the random attenuation and <^(r) is the random phase noi.se process. If 
ower fading phenomena (power law attenuation, and log-normal fading (Rappaport 
) are compensated for by power control, and the multipath phenomenon is spatially 
)genous, then the process rit) is stationary with a Rayleigh marginal distribution R, 
density 


PRir) = -jexp 





We note that E{R}) = 2a^. The time variations in the process r{t) are of the order of the 
Doppler frequency fd, which is related to the carrier frequency /c, and the vehicle speed 
V, by the formula 

fd = —, ( 2 ) 

c 

where c is the speed of light. For a 900 Mhz carrier frequency, for example, the above 
formula yields a Doppler frequency of 3 Hz/m/sec, or 10.8 Hz/km/hr. Thus, for the sig¬ 
nalling rates used in high speed wireless transmission (e.g.. Mbps), Rayleigh fading can 
be taken to be roughly constant over several bits (see also Rappaport 1996, pp 165-166). 

From (1), the predetection signal-to-noise ratio (SNR), say ifr (r), at the receiver is given 
by 


fit) = (r(0)^ 



xrnit 


( 3 ) 


where (Eh/ !^Q)xmit is the “transmitted” SNR. Denote the marginal random variable for 
the stationary process f(t) by 4'; then we have 

^ = R^iE,,/No)„ni,, (4) 


where R is the marginal of the process r{t), as defined above. The average effective 
received SNR is then given by 



In dB units we can now write (4) as 


(4')i/B = r^) +iA)dB (5) 

\NoJdB 

where A := R^jla^, and the distribution of A is known to be given by P(A > a) = 
(Rappaport 1996). As observed above, the SNR f(t) varies slowly as compared to the 
signalling rate. When the SNR is low (i.e., a large negative values of iA)dB)> this situation 
persists for a while and the bit error rate (BER) is high; then the SNR improves and stays 
high for a while, yielding a low BER. This view motivates a Markov model for packet 
loss; we develop the model in the next subsection. 


2.2 The Markov packet loss model 


For a given average SNR (Eb/Nq) (say 30 dB) we determine the amount by which the 
SNR must drop below the average so that the channel enters the “Bad” state (say 20 dB, 
i.e., a factor of .01); call this “margin” a, in dB units. Then, defining 5 := 10"“'^’*^ (i.e., 
S = .01 for a = 20 dB), the fraction of time that the channel is in the Bad state is given by 
P(A < 5) = 1 - e~^. Thus, fixing a fixes the fraction of time that the channel is in the 
Bad state. To obtain the mean durations in each state, we use results from level crossing 
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(1-P) 


Figure 2. Transition structure of the Markov 
loss model. 


/sis of the process r{t) (see, Parsons 1992). Defining G(5), and 8 ( 8 ), as the mean 
tions in the Good and Bad states, respectively, the following formulas are obtained 


G(S) = 


8 ( 8 ) = 


fd 


-1 


•v/2i5’ 


1 ) 


e fd is given by (2). For small S, 
8 ( 8 ) ^ 8 G( 8 ). 


1 Rs S, hence (6) and (7) yield 


( 6 ) 


(7) 


( 8 ) 


loss model is a discrete “time” Markov chain whose transitions are embedded at packet 
idaries; thus we express the Good and Bad periods in units of packet transmission 
on the link. The transition structure of the model is depicted in figure 2. The chain 
sumed to be embedded at the beginnings of packet transmissions. Thus, if a number 
ickets are being transmitted back to back, and if the chain is in the Good state when 
:ket is about to be transmitted then this packet will be Good, and the next packet will 
ood or Bad with probabilities 1 — y and y, respectively. It follows, from (8), that 



ir calculations we will use the above approximation, as S <0.1 (i.e., a > 10 dB) in 
mmerical results. Observe that we are making the tacit assumption that the durations 
ig which the SNR is above or below the threshold level ai'e exponentially distributed; 
Is not true (see Rice 1958), but the minimal Markov model that we have considered 
ot model any additional characteristic of the fading process. 

Lven a, equivalently 8 , the individual values of y and 5 in the Markov model are 
ined from the Doppler frequency fd (2), and (6) and (7). For example, consider packet 
mission times of 7.5 ms (e.g., 1.5 Mbps link, and 1400byte packets), fc = 900MHz, 
a = 20dB. For a speed of 5kmph, fd — 4.17 Hz, the mean good period is about 
packets, the mean bad period is about 1.28 packets, and y = 0.0078; at a speed of 
cmph, fd = 83.3 Hz, the mean good period is about 6.4 packets, the mean bad period 
out 0.064 packets, and y = 0.156. Thus, clearly, some care is needed in using this 
el for high vehicle speeds at which the fade duration is smaller than one packet length. 
lOut a link protocol, whenever there is a fade during a packet transmission, that packet 
5t, even though the channel is good during some parts of the packet. For high speeds, 
/hich the fade durations are smaller than a packet (but not so high that G + 8 < mean 


1 iO 


/\nuru^ i\uinur unci juck, nuiii.rnun 


pkt xmission time), we take = 1 (i.e., exactly one packet is lost in each fade), and y is 
obtained from the following equation 

G + B _ 

mean pkt xmission time 

Thus, at low speeds, the probability that a packet is bad is just 5/1 +- S, whereas for high 
speeds, when the fades are shorter than the packet transmission time, the probability that 
a packet is bad increases with increasing speed, even though 8 is fixed. 




2.3 The loss model with a link protocol 


The round trip propagation delay on teiTestrial mobile radio links is typically smaller than 
the packet transmission time. Consequently, a stop-and-wait link protocol suffices. We 
characterise a link protocol by the link packet length and the number of times it attempts 
a packet. 

Observe that with 128byte link packets and 1.5 Mbps transmission, the link packet 
transmission time is about 0.68 ms. For 1400byte TCP packets, there are 11 link packets 
in a TCP packet and for speeds above 30 or 40kmph the mean fade duration is less than 
2 or 3 link packets. Thus, with a link protocol, we embed the Markov packet loss model 
(figure 2) at the beginnings of link packet transmissions. We use the loss model to obtain 
the probability of TCP packet loss, and assume that TCP packet losses are independent 
at the TCP packet level. This assumption is adequate if the link packets are small and the 
number of link packets in a TCP packet is large. We assume that the first link packet in a 
TCP packet finds the Markov model in its stationary distribution. The probability of TCP 
packet loss is then the probability that for at least one link packet all attempts fail. The link 
protocol also causes an increase in TCP packet transmission time; we use the link level 
Markov loss model to obtain the inflated mean TCP packet transmission time. 

Let N denote the maximum number of attempts of a link packet; and L denote the 
number of link packets in a TCP packet. Recalling that the Markov process is embedded 
at the beginnings of link packet transmissions, we define 


Pn = Prob{at least 1 out of n link packets fails, given that the initial state is Good} 

<?n*^ = Prob(at least 1 out of n link packets fails, given that the first link packet has already 
had k(< N) bad attempts and the initial state is Bad), 

p = Prob{a TCP packet is bad). 

Assuming that the first link packet in a TCP packet finds the Markov process in its 
stationary distribution, we have 


P = 




y + ^ 


PL + 


y 


y + P 


d- 


The values of, 1 < n < L and gn^\ l<n<L,0<k<N—l are easily obtained from 
the following equations by recursive substitution, with the boundary conditions: pi = 0, 


(N- 

Qn 


1) 


A <n < N, and qf'^ = (1 — ^\0 < k < (N — 1). For n > 1 and 
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:k < N 

Pn = (i -y)p„-i + y<7^-i> 
qjP = fip„ + (\ - 

obtaining the mean effective TCP packet transmission time, taking into account the 
c retransmissions, define 


„ = mean number of link packets that will be sent for transmitting a TCP packet of 
length n link packets, given that the state at the beginning is Good, 

^ = mean number of link packets that will be sent for transmitting a TCP packet of 
length n link packets, given that k link packets have already been sent for the first 
link packet and the cuiTent state is Bad, 

/ = mean number of link packets that will be sent for transmitting a TCP packet of 
length L link packets. 

As for the loss probability, we have 


/ = 




Il + 


y 

y 


m 


0 

L- 


;h the boundary conditions /i = 1 and m\^ = 1, we have the following equations, 

n > 1 and 0 < k < N — I, 

= 1 + (1 + 

4, = 1 + (1 — y)ln-\ + y^n^^l.\-, 

jse equations can also be solved recursively. 


TCP window adaptation 

any time t there is a lower window edge A{t), which means that all data numbered 
to, and including A(t) — I has been transmitted and acknowledged. The transmitter 
send data numbered A(t) onwards. There is also the transmitter’s congestion window 
t), which, at time t, is the maximum amount of data that the transmitter is permitted to 
d, starting from A (r). Under normal data transfer, A(0 has nondecreasing sample paths, 
ereas the adaptive window mechanism causes W (t) to increase or decrease, but never 
eed Wmax- The receipt of an acknowledgement (ack) packet that acknowledges some 
a will cause a jump in A(0 equal to the amount of data acknowledged, but the change 
W (t) would depend on the particular version of TCP and the state of the congestion 
Ltrol process. 

U the transmitter, round-trip time measurements of some packets are used to obtain 
Linning estimate of the packet round-trip time (rtt) on the connection. Each time a 
v packet is transmitted, the transmitter starts a timer and resets the already running 
ismission timer, if any. The timer is set for a round-trip time-out (rto) value that is 


derived from the rtt estimation procedure. Whenever a time-out occurs, retransmission is 
initiated from the next packet after the last acknowledged packet (i.e., from the sequence 
number Ait)). It is important to note that the TCP transmitter process measures time and 
sets timeouts only in multiples of a timer granularity, for example, BSD UNIX based 
systems have a timer granularity of 500 ms. Further, there is a minimum timeout duration 
of 2 or 3 timer “ticks” in most implementations. We will see, in the analysis, that coarse 
timers have a significant impact on TCP performance. For details on rtt estimation, and 
the setting of rto values, see Desimone (1993) or Stevens (1994). 

The following basic window adaptation procedure (Van Jacobson 1988) is common to 
all the TCP versions; our description follows that of (Lakshman & Madhow 1997). At all 
times t the following are defined for all the protocol versions: 

W(0 = the ti'ansmitter’s congestion window at time t, 

Wihit) = the slow-start threshold at time t. 

The evolution of these processes is triggered by acks (only first acks; not duplicate acks) 
and timeouts as follows. 

1. If W(t) < Wth (t), each ack from the receiver causes W(t) to be incremented by 1. This 
is called the slow start phase. 

2. IflV(f) > VP(/,(r),eachackfromthereceivercauses VP(0tobeincrementedby l/lT(r). 
This is called the congestion avoidance phase. 

3. If timeout occurs at the transmitter, at epoch t, Wif^) is set to 1, Wthif^) is set to 
rW(f)/21, and the transmitter begins retransmission from the packet after the last 
packet acknowledged. 

Note that the transmissions after a timeout always start with the first lost packet. We 
will call the “window” of packets transmitted from the lost packet onwards, but before 
retransmission starts as the loss window. 

If a packet is lost, then the ack number from the receiver will cease to be advanced, and 
the transmitter at the LAN host will continue sending packets until the current window is 
exhausted, or a packet gets through to the receiver and an ack is received. The last packet 
transmitted will have a timeout associated with it; in TCP-OldTahoe, retransmission will 
start only upon the expiry of this timer. Later versions of TCP, i.e., Tahoe, Reno, and 
NewReno, implement a fast-retransmit scheme. Observe that, even if a packet is lost, 
if subsequent packets get through, the receiver will continue to send back acks, but the 
“expected packet” number is not incremented. If several (an integer parameter K, e.g., 3) 
such “duplicate” acks are received at the LAN host, then it can assume that the loss is 
a random loss, and not due to congestion. When the number of duplicate acks received 
reaches the threshold K, then the packet next expected by the receiver is retransmitted. A 
complicated recovery phase then follows. In particular, in TCP Tahoe (Fall & Floyd 1996) 
if the transmitter receives the K th duplicate ack at time t, before the timer expires, then the 
transmitter behaves as if a timeout has occurred and begins retransmission, with W(r”*") 
and W//j(r+) as given by the basic algorithm. 
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. TCP throughput analysis with the Markov loss model 

.ssume that, at t = 0 the connection starts with W(0) = 1 and W,i,{0) = rM^max/21, 
here Wmax is the maximum window limit set by the receiver; this usually depends on 
le receiver’s buffer size (typical values are 32 kbytes or 64 kbytes, i.e., several lO’s of 
CP packets). The protocol starts at time 0 in slow start. Let be the epoch at which 
le first packet loss occurs, and let U] = W(£i). As described above, timeout or fast 
jcovery follows and at some later epoch transmission resumes with the first lost packet, 
’this epoch finds the wireless channel in its Bad state, then a loss occurs immediately, the 
meout is doubled, and the next cycle starts with W = 1 and Wth = 2, the minimum value 
f Wth (Stevens 1994). No successful packet is sent in such a period. Hence, we merge 
ich periods, if any, into the recovery period of the previous cycle. Thus at the beginning 
f the next cycle, denoted by t \, the channel is in the Good state. Again an epoch ^2 can be 
lentified at which the first loss occurs in this transmission cycle. For ^ > 1, let 4 denote 
le kth epoch at which a new transmission cycle starts as described above. For k > I, 
'e call the interval ( 4 _i, tk] the ^th cycle. In the A:th cycle, let £k be the epoch at which 
le first packet is lost in the cycle (this is an end of a packet transmission epoch from the 
)uter). Further, fork > 1 , let = W {lu) denote the transmitter’s window size at which 
acket loss takes place. We take Uq — VFmax- Since the first packet transmitted in each 
/cle is always good, the state space of the process {fy^} is (2, 3,..., VKmax}- 

If the first retransmission after a recovery period finds the channel in the good state, 
len for TCP-OldTahoe and Tahoe, the value of Uk determines the values of W(t^) and 
according to W(t^) = 1 and Wthit'^) = rt/A:/21. We know that at the beginning 
f the transmission of the packet that is lost, the channel is in the Bad state. It is thus clear 
lat the evolution of the congestion window process after the first loss epoch in the kth cycle 
spends only on the value of f/*; thus the process {Uk} is a discrete time Markov chain 
ver the state space {2,..., Wmax}- We can obtain the stationary distribution of this chain, 
fote that a more complex analysis ensues if losses can occur in either state of the channel. 

The mean number of packets successfully transmitted in each cycle, before the first lost 
acket, is just 1/y. When a loss occurs at the lossy link, there would be some packets 
aeued in the lossy link transmitter and the TCP transmitter at the LAN host will continue 
mding packets till the congestion window (Uk) is exhausted (even if fast-retransmit is 
nplemented, by the time 3 duplicate acks are received, a fast sender would already have 
^hausted the window). Some of the packets that are transmitted on the lossy link, after 
le first lost packet, will get through, and will not need to be resent. We denote the mean 
amber of such packets as residuaLgood. 

Thus, for each TCP version, we have a Markov renewal reward process (Wolff 1990), 
aibedded at the epochs tO) fh fz, • • •, the reward being the number of good packets put 
irough in each cycle. Then the throughput T is given by 

y, ^ 1/y + residuaLgood 
mean_cycle-duration 

With a link.protocol and high vehicle speeds, we assume a Bernoulli loss model at 
le TCP packet level. Hence, for this case the analysis in (Kumar 1996) carries through 
irectly, after accounting for the fact that TCP packet transmission times are longer. 


We will assume that the minimum timeout is large compared to the other time scales 
in the local network; this is true for the numerical parameters that we used (Kumar 1996) 
(0.75 s minimum timeout, as in BSD and 7.5 ms packet transmission time). Then, if the 
loss in a cycle is followed by a timeout, we assume that the beginning of the next cycle 
finds the Markov loss process in its stationary distribution. This is always the case for 
OldTahoe, since it recovers only after a timeout. In the case of Tahoe, however, if fast 
retransmission takes place then we need to know the state of the Markov loss process at 
the beginning of the next cycle. 

The following analysis also assumes that the LAN host can transmit much faster than 
the wireless link. Hence, during the packet transmission periods in a cycle, the wireless 
link is assumed to be continuously busy. This assumption facilitates the carrying over of 
the channel state from one packet to the next. 

4.1 Analysis of the markov chain [ Uk ] 


Owing to the above viewpoint, the transition probabilities of {Uk} are slightly different 
from the ones in our earlier paper (Kumar 1996). For M G {2,..., Wmax). define 

9m = Prob{the next cycle starts with Wth = \M/T\, given that the loss window in this 
cycle is M] 

9m = Prob{the next cycle is forced to start with Wth = 2, owing to additional timeouts, 
given that the loss window in this cycle is M (= 1 — 0 a/) }. 


Following the earlier development (Kumar 1996) we can now write down the transition 
probabilities. Given that f/^ = M G {2,..., Wmax} (and m := max(2, [M/2])), recalling 
the Markov loss model, and writing y = 1 — y, define two transition probability matrices 
P and Q as follows: 


PMJ = 


yO‘-l)-ly 


P3,ji= P2J=PAj) 


2 < j < m — ] 
j = m + k, 

^ S — (l^max — 1 — m) 

J = Wmax 

for 2 < j < Wmax 


where 


^(fc) = (fc+l)((m-l) + (fc/2,..)). 

Note that P con*esponds to the case in which the next cycle starts with Wth = \Ml/T\ and 
Q corresponds to the case in which the next cycle is forced to start with Wth = 2. Finally, 
we define the vector 


0 = (( 92 , 

and then the transition probability matrix of {Uk} is written as 

diag(0)P + (/- diag(0))!2> ( 10 ) 

where / is the (Wmax — 1) x (Wmax — 1) identity matrix. 



Performance analysis of TCP over a mobile radio link 


123 


The mati’ices P and Q are common to the OldTahoe and Tahoe versions with the same 
lines of ITmax and Markov loss model transition probabilities y and f ; as shown below 
depends on y and f for both OldTahoe and Tahoe. 

IdTahoe: Since we assume that after a timeout the Markov channel model is found in its 
ationary distribution, for M € {2,..., M^max}. 


9 m = 


Y + fi 


= 1 — 6m- 


ihoe: Define, for M € {2,..., H/max). 

bM = Prob{the good transmission period in a cycle is followed by fast retransmission, 
if the loss window is M] 

' r'\ 

— Prob{the good transmission period in a cycle is followed by fast transmission, and 
at this epoch the channel is in the Good state} 

= Prob{the good transmission period in a cycle is followed by fast transmission, and 
at this epoch the channel is in the Bad state (i.e., = <pM. “ 

follows that 

+ (1 — 

Recursive equations were developed for computing and (Pm '' owing to lack of 
lace here we refer the reader to the full report (Kumar & Holtzman 1996). 

Thus the transition probabilities of the process [Uk] can be calculated by first com- 
iting (pM< (P)w^<4>m \ and 9m, and finally using (10). The stationary distribution can 
; obtained by any of the many standard techniques. Denote the stationary distribution 


2 < M < Wj-nax' 

2 Throughput analysis 

iven that the loss window is M, the probability that k < {M — 1) of the remaining packets 
:t through isjustg^^_, + b^M_y Given the stationary distribution ha/, 2 < M < Wmax, 
e mean number of good packets transmitted in a cycle after the first loss occurs can thus 
! calculated. This is what we had called residuaLgood (see (9)). 

We adopt the approximate mean cycle time analysis discussed in § 4.2 of Kumar 
996). Letting Z denote the mean of the loss window distribution, we take the aver¬ 
se window during packet transmissions to be 0.75Z, and take the network throughput 
be r = r(A, 0.75Z) (A is the LAN packet transmission rate normalised to the wire- 
ss link transmission rate), where r(X, w) ~ (7.^"^+^^ — A)/(A.^'"‘'‘’' — 1), if A. 1, and 
1, w) = u)/(l -1- w). 

We assume that the timeout is always the minimum timeout (rfo-min) (a good as- 
imption for a large coarse timeout and the local network situation). It follows that the 


recovery time due to repeated, exponentially growing timeouts, is given by rto_ min/ 
[1 — 'lyjiY + ^)] (for this equation to hold, we must have y/{y + yS) < 0.5). 

It follows that the throughput of OldTahoe is given by 

^ _ (1/y) + residual-good 

^ (l/yr) + rfo-min/[l-2y/(y + j6)]’ 

where the first term in the denominator is the mean time for transmitting 1 /y good packets. 
The throughput of Tahoe is given by 


T’oidTahoe = (1/y + residuaLgood) 


M=2 


rto-mm 


2(y/(y + /3)) 


^(G) 


M 


+ <pMy + ^M 


(fi) / M 


7 + 


rto ^rmn 


l-2y/(y4-)S), 


Here the second tenn in the denominator is an expectation over the loss window distri¬ 
bution; the term conesponding to each M is the sum of the recovery durations for three 
possibilities: the first, if fast retransmit fails, the second, if fast retransmit succeeds and 
the next cycle finds the channel in the Good state, and the last, if fast retransmit succeeds 
but the channel is found in the Bad state, necessitating a timeout-based recovery. 


5. Numerical results and discussion 

We present numerical results for the same numerical parameters as in Kumar (1996). The 
wireless link speed is 1.5 Mbps (all rates are normalised to this rate) and the LAN host can 
transmit at 5 times the wireless link rate (i.e., L = 5). The TCP packet length is 1400 bytes; 
i.e., its transmission time on the wireless link is 7.5 ms;.various times are normalised to 
this transmission time. The minimum timeout is 750 ms; or 100 packet service times on 
the wireless link. The fast retransmit threshold is K' = 3. The carrier frequency fc is taken 
to be 900 Mhz. 

We consider DPSK (differential phase shift keying) modulation. In this analysis we do 
not consider forward error correction coding, or bit interleaving. The link protocol operates 
directly over the modulation scheme, and has an error detection code. With AWGN the 
BER for DPSK is given by 

f = 0.5exp[-(£i,/iVo)]. (11) 

Thus without a link protocol, in AWGN, the TCP packets can be assumed to be lost 
independently with probability p given by 

p = 1 - (1 - e)PacketJength ^j2) 

We will provide numerical results with AWGN for comparison purposes. In figure 3 we 
plot the throughput of TCP OldTahoe and Tahoe versus the average SNR, with and without 
fading (i.e., AWGN). There is no link protocol; thus, whenever a packet encounters a fade 
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TCP Thruput vs. Mean SNR, Wmax=24, Iink1= 5 x Iink2, min t.o.=100 



Figure 3. Throughput of versions of TCP vs. mean signal to noise ratio; no link protocol; 
K is the fast retransmit threshold. 


at packet is certainly lost. In the case of fading, two curves are shown for each version, 
responding to average fade durations of 1 and 2 packets. The Bad state corresponds to 
iSNR< 10 dB. 

Figure 3 shows that with AWGN an SNR of about 12dB yields very good throughput, 
hereas, with fading, an SNR of about 30 dB is required to obtain a throughput better 
an 90% from TCP Tahoe. TCP OldTahoe requires an SNR of 40dB. For a given SNR, 
creasing the fade length appears to improve the TCP throughput. This observation is 
iplained as follows. The “fade limit” of 10 dB, taken together with the average SNR, 
tes a (see Section 2.2), and hence fixes (5; i.e., for a given value of Ei,/No the ratio of 
e good and the bad periods is fixed (see (8)). It follows that, for fixed Eh/No, if the 
de duration is increased from 1 to 2 packets then the good periods also increase by a 
ctor of two. Thus, although increasing the fade duration results in a greater frequency 
' consecutive losses, since the good duration also increases, the TCP window can build 
large values before losses do occur. This yields a larger throughput for increasing 
de duration. These observations hold for low speeds at which the fade durations are 
)mparable to the packet transmission times (e.g., E^/Nq = 20 dB implies speeds of 
)out 10 to 20kraph). 

Consider what would happen if the vehicle speed was allowed to reduce to zero. For 
:ample, if Ei,/Nq = 30 dB then the probability that a connection starts in the bad state 
0.01 (see (8)). For very low speeds, during the entire duration of the connection, the 
lannel will be in the same state as it was found at the beginning of the connection. Thus 
the initial state found is bad the throughput is 0, otherwise the throughput is 1. Hence, 
'eraged over connections, the average throughput seen by a TCP connection will be 
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0.99. It follows that in figure 3, for a fixed SNR of 30 dB, as the vehicle speed decreases 
the throughput of both Tahoe and OldTahoe will converge to 0.99. Similar calculations 
can be done for each value of SNR; these results are plotted as “limit, speed 0” in 
figure 3. Note that we have no link protocol in the model for which results are shown in 
figure 3. With a link protocol, the throughput should be no worse than without one. Thus, 
we see that for speeds approaching zero the average throughput increases and converges to 
the throughput obtained from a “quasistatic” analysis assuming veiy long fade durations 
and good durations but with the appropriate probabilities. This observation is consistent 
with the results with a link protocol that we present next. We now fix the average SNR 
(to 30 dB or 25 dB), and, for a fixed fade limit of 10 dB, we study TCP throughput as a 
function of vehicle speed. We now include the link protocol model in our analysis (see 
Section 2.3). 

In figures 4 and 5 we show the TCP throughout versus vehicle speed. The link packet 
length is 128 bytes; thus there ai*e 11 link packets in a TCP packet. In figure 4 each link 
packet is attempted N = 2 times, whereas in figure 5 each packet is attempted N = 8 times. 

Figure 4 shows that, as the vehicle speed increases from a small value, first the throughput 
(for both versions) decreases, hits a minimum, and then it increases. This is explained as 
follows. Note that if the fade duration is longer than the number of attempts made for a link 
packet then that link packet and the corresponding TCP packet is surely lost. To the left 
of the minimum point, the mean fade duration is longer than 2 link packets, and increases 
as the speed decreases; however, as observed earlier, the good period length increase too, 
hence the throughput increases. To the right of the minimum point, the fade duration is 


TCP Thruput vs. Vehicle Speed, Wmax=24, Iink1 = 5 x linka, min t.o.=100 



Figure 4. Throughput of versions of TCP vs. vehicle speed, for fixed Eh/N q = 30 dB and 
fade limit = 10 dB; each link packet is attempted 2 times. 
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TCP Thruput vs. Vehicle Speed, Wmax=24, Iink1= 5 x Iink2, min t.o.=100 



Figure 5. Throughput of versions of TCP vs. vehicle speed, for fixed Eb/N{] = 30 dB, 
and fade limit = 10 dB; each link packet is attempted 8 times. 


TCP Thruput vs. Vehicle Speed, Wmax=24, Ilnk1= 5 x ljnl<2, min t.o.=100 



TCP pkt=1400B, link speed=1.5Mbps, llnkpkt=:128B, 8 attempts 

Figure 6. Throughput of versions of TCP vs. vehicle speed, for fixed £/,/A^o = 25 dB, 
and fade limit = 10 dB; each link packet is attempted 8 times. 
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smaller than 2 packets, and for larger vehicle speeds there is increasing probability that 
not all attempts of a link packet fail. 

We note that the behaviour of throughput versus vehicle speed, for fixed SNR, in 
figure 4 is consistent with our discussion for low speeds in the context of figure 3. The 
throughput does increase for low vehicle speeds. If the curves in figure 4 are extrapo- 
lated as speed goes to zero, however, the limit will not match that obtained in the con¬ 
text of figure 3. This is because the model used to obtain the results in figure 4 (fading 
model at the link packet level, but Bernoulli losses across TCP packets) does not apply 
for very low vehicle speeds for which there is significant fade correlation between TCP 
packets. 

Figure 5 shows that each link packet needs to be attempted 8 times for Tahoe to yield 
a good throughput for a large range of vehicle speeds. In figure 6 we show the effect 
of reducing the SNR by 5dB, We find that 8 attempts are no longer enough even for 
Tahoe. The number of link packet attempts need to be increased to 16 for Tahoe to provide 
reasonably good performance (see Kumar & Holtzman 1996). 


6. Conclusion 

For the default parameters of the BSD implementation of TCP, over a 1.5 Mbps wire¬ 
less link, as a general observation, we find that (without physical level enhancements, 
such as forward error correction and bit interleaving*) an SNR of 25 dB to 30 dB is re¬ 
quired to obtain from Tahoe a throughput greater than 90% of the wireless link speed. 
Without a link protocol, OldTahoe requires 40 dB; and with a link protocol OldTahoe 
suffers at low vehicle speeds, for which the fade durations are large. When using a link 
protocol, the choice of link packet length and the number of attempts for each packet 
are important parameters. Basically, the attempts have to “outlast” the fade. Since very 
large link packets defeat the idea of using link packets, small link packets have to be at¬ 
tempted several times. For a fixed SNR, a given link packet length, and a given number 
of link packet attempts, we find that throughput varies in an interesting way with vehicle 
speed. At very low speeds (e.g., pedestrian speeds) the throughput is high; it decreases 
with increasing vehicular speed until the fade durations become shorter than the num¬ 
ber of link packet attempts. Beyond this speed, again the throughput increases, and is 
limited only by the expansion of the TCP packet transmission time owing to link-level 
retransmissions. 

Much work remains to be done to develop more comprehensive protocol analyses even 
with the simple two state Markov loss model. Our analysis at present does not apply 
to many situations of interest; e.g., link protocols at very low speeds, or with large link 
packets. The Markov model can also be enhanced to include a nonzero loss probability 
in the Good state. Both these situations will require a more complex analysis, as the state 
of the Markov loss model will need to be included in the evolution of the TCP window 
process. Such enhancements of the analysis will be fruitful as future research. Other TCP 
versions, such as Reno and NewReno, also need to be analysed. 


In any case, forward error correction and interleaving are not much help when there is slow fading. 
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